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1 . INTRODUCTION 

For  the  last  two  years  A.R.A.P.  has  been  developing  a 
computer  program  capable  of  solving  the  three-dimensional 
steady  flow  problem  of  fluid  motion  in  a stratified  fluid. 
This  work:  was  accomplished  in  various  stages,  and  is  now 
assembled  in  one  FORTRAN  program  called  WAKE.  This  program 
resides  on  the  A.R.A.P.  computer  system.  This  part  of  the 
final  report  summarizes  the  usage  of  WAKE  and  gives  detailed 
explanations  of  its  input  requirements  and  output  results. 

In  Section  2 we  briefly  summarize  the  theoretical  problem. 

In  Section  3 we  examine  the  numerical  scheme  used  to  solve 
the  equations  of  motion.  In  Section  4 we  detail  the  input/ 
output  specifications,  including  the  structure  of  the  impor- 
tant initial  profile  file.  The  Appendix  gives  a source 
listing  of  the  entire  .JAKE  program  and  its  full  subroutine 
complement.  This  guide  is  not  intended  as  a full  explana- 
tion of  the  WAKE  program  - only  the'  FORTRAN  listings  can  do 
that.  Nor  Is  it  intended  to  demonstrate  the  conversion 
possibilities  of  WAKE  to  other  computer  facilities.  Rather*, 
with  this  guide  a computer*  analyst  unfamiliar  with  the  WAKE 
program  should  be  able  to  construct  a needed  set  of  initial 
profiles,  grasp  an  overview  of  the  numerical  scheme  and 
program  structure,  and  produce  suitably  correct  output  from 
the  A.R.A.P.  computer  facility  in  a fairly  short  period  of 
t i me . 
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2.  THE  THEORETICAL  PROBLEM 


When  a sell‘-propelled  body  moves  through  a medium  with 
a stratified  density  gradient,  a wake  is  generated  which 
expands  behind  the  body  as  its  potential  energy  increases. 

Ear  behind  the  body  the  potential  and  kinetic  energies  come 
into  balance  and  the  wake  collapses.  The  heavier’  fluid, 
finding  itself  in  a region  of  lighter  background  density, 
reestablishes  the  stable  condition  existing  before  passage 
of  the  body,  but  at  the  expense  of  the  generation  and  trans- 
mittal of  internal  gravity  waves.  The  theoretical  prediction 
of  this  complex  physical  problem  has  been  the  subject  of  a 
great  deal  of  study.  The  intent  of  the  A.R.A.P.  approach  is 
to  model  the  generated  turbulence  by  the  technique  of  invariant 
second -order  closure  and  to  follow  the  buildup  and  collapse 
phase  through  two  Brunt -Valsal'a  (B.V.)  periods  of  fluid  motion. 
A more  detailed  explanation  of  the  derivation  of  the  equations 
Is  presented  in  Part  1 of  this  final  report  (ref.  1).  For 
completeness  we  here  present  the  derived,  modeled,  approximated 
and  normalised  equations  as  they  stand  prior  to  numerical  solu- 
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For  the  perturbation  density  p: 
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For  the  mean  velocities  u,  v,  and  w; 
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For  the  scale  length  A: 
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For  perturbation  pressure: 
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For  the  turbulent  correlations: 
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A complete  nomenclature  may  be  found  in  Part  1.  Restric- 
tions and  modifications,  especially  in  regard  to  the  Quasi - 
Equilibrium  eqs.  (2.8)  for  the  turbulent  correlations  and  the 
correction  factor  f,  are  also  detailed  in  Part  1.  The  algebraic 
solution  to  eqs.  (2.8),  with  the  assumptions  that  the  principle 
production  gradients  are  in  u and  and  density  gradients  in  v 
are  smaller  than  density  gx-adients  in  z (for  computation  of  f 
only)  gives: 
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The  solution  procedure  always  involves  the  solution  ol'  eqs.  (2.9) 
for  the  turbulent  correlations  and  eq.  (2.6)  for  the  dynamic  scale  A. 

o 's 

However,  the  solutions  of  the  main  variables  q , , , u,  v,  w,  and  7/ 
are  governed  by  the  regime  of  interest.  Thus,  only  q‘ , p and  u are 
computed  in  a Phase  I calculation  where  R1q  <0.1.  Here  v = w = v - 0 

since  the  flow  is  far  from  collapse.  In  a Phase  II  calculation,  we 

2 

solve  for  q , p,  v,  w,  and  v,  assuming  that  u = 0;  we  are  here 
restricted  to  flows  for  which  |u  I /q  <0.1.  When  the  flow 

situation  does  not  fit  either  condition,  we  calculate  all  the 

2 A 

variables  q , p,  u,  v,  w,  and  u in  a Phase  III  calculation. 

For  the  runs  presented  in  Part  I (ref.  l),  we  take  our  initial 
conditions  on  the  main  variables  (when  they  are  nonzero)  to  be  those 
assembled  in  Table  2.1. 
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TABLE  2.1 

TABULATION  OP  NONZERO  INITIAL  CONDITIONS 
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3.  THE  NUMERICAL  APPROACH 


Equations  (2.1)-(2.8)  are  solved  in  the  two-dimensional 
y -z  plane  by  recasting  them  in  finite  difference  form  and 
applying  the  ADI  technique  of  Peaceman,  Rachford  and  Douglas 
(refs.  2 and  3).  In  this  method  we  construct  a two-dimensional 
grid  in  the  y and  z directions  and  march  in  the  x direction  as 
we  follow  the  flow  development  downstream  of  the  initial  condi- 
tions. Tn  the  y and  z directions  the  first  and  second  deriva- 
tives are  approximated  by  centered  differences  (spacing  is 
variable),  while  forward  differences  are  used  in  x (ref.  4). 
Inputted  tolerance  parameters  control  the  size  and  speed  of 
the  marching  direction,  and  the  spacing  variability  and 
intensity . 


At  the  beginning  of  a new  step  in  Ax  (perhaps  the  start  of 
the  run  itself),  the  program  performs  an  implicit  sweep  in  one 
direction  in  Ax/2  and  then  sweeps  in  the  other  direction  in 
Ax/2  to  complete  the  ADI  procedure.  The  initial  sweeping 
direction  alternates  with  each  full  step  to  unblas  any  solution 
near  the  edge  of  an  expanding  profile.  These  steps  have  used 
the  current  values  of  q‘ , , , u,  v,  w,  and  h , together  with  the 
gross  scale  A to  step  forward  in  x.  Although  the  main  variables 
are  coupled  by  eqs.  (2.1)-(2.5)  and  (2.8),  we  choose  to  use 
current  values  wherever  necessary  to  decouple  the  equations 
completely . solutions  at  the  next  x are  then  swept  to  compute 
maximum  values,  maximum  changes,  and  various  integrals  of 
interest,  along  with  the  updated  scales  A and  \ . The  next 

v Z 


step  Ax^ew  is  computed  based  on  the  changes  taken  by  the  present 


step  in  relation  to  the  maximum  change  permitted.  The  profiles 
are  then  swept  again  to  update  the  pressure  forcing  function  - 
the  right  side  of  eq.  (2.7) > and  to  compute  the  algebraic 
turuulence  via  eq.  (2.8). 
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WAKE  then  calls  the  pressure  Iteration  subroutine.  It  works 
as  a miniature  WAKF  program  by  adding  a - ch/dt  term  to  the  left 
side  of  eq.  (2.7)  and  performing  iterations  in  y stepping  from 
the  current  solution  to  the  next  steady  state  estimate.  Appro- 
priate output  routines  are  then  called,  followed  by  the 
necessary  set  of  routines  that  inspect  the  curvatures  in  the  y 
and  z directions,  readjusting  the  profile  again  per  inputted 
tolerance  criteria.  A new  step  is  then  taken. 

The  mainline  program  is  called  WAKE.  The  Appendix  of  this 
report  contains  a complete  listing  of  WAKE  and  its  subroutines; 
in-house  disk  I/O  routines  and  other  straightforward  assembler 
routines  are  not  included.  These  routines  control  the  monitoring 
of  the  disk  files  storing  the  large  array  of  data  necessary  to 
execute  the  calculation. 

It  may  be  worthwhile  to  realize  that  the  file  record  length  is 
24  words,  broken  into  four  sections  of  six  words  each.  The  first 
section  contains  TV(o),  where  TV(l)  ^ du/dy;  TV(2)  = du/dz; 

TV(3)  = dp/dy;  TV(4)  = dp/dz ; TV(t>)  = F;  and  TV(b)  ^ y . The  second 
section  holds  XV(6)  at  the  present  step  value;  the  third  holds  the 
intermediate  (first-sweep)  values  of  XV(6);  and  the  fourth  section 
holds  the  new  step  values  of  XV(6).  Here  XV(l)  - q2,  XV(2)  = p; 
XV(3)  = u;  XV( 4 ) - v;  XV(b)  w;  and  XV(o)  - A. 
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4.  INPUT  CONTROL 


Initialization  of  a WAKE  ran  requires  the  generation  of 
an  initial  profile  file  (giving  the  desired  initial  conditions 
to  the  desired  variables)  and  the  input  of  an  appropriately 
punched  deck  of  computer  cards.  The  initial  file  must  be  formed 
in  a way  totally  compatible  with  the  sample  generation  program 
shown  at  the  end  of  the  Appendix.  It  must  be  structured  as 
follows . 


First  word: 

Second  word: 

Next  NY  words  : 

Next  NZ  words : 

Next  NY*NZ  words: 


Next  NY*NZ  words: 


the  initial  x position  value. 

the  number  of  points  NY  in  the  y direction  and 
NZ  in  the  z direction. 

the  NY  independent  y values  of  the  y direction 

mesh  (monotonically  incre  sing). 

the  NZ  independent  z values  of  the  z direction 

mesh  (monotonically  increasing). 

the  complete  corresponding  mesh  values  for 

the  first  Initialized  dependent  variable. 

The  file  must  contain  these  NY*NZ  values  in 
blocks  of  10  y values  at  a time  (for  all  z), 
until  the  last  block  contains  enough  values 
to  reach  NY.  Thus,  the  blocks  would  be  built 
as  10*NZ , 10*NZ , 3*NZ  if  NY=25. 
the  second  initialized  variable. 


Next  word:  - 1.0  to  signal  end-of-data. 

With  the  disk  buffering  currently  in  operation,  the  file  inver- 
sion program  PBFFI  must  be  called  to  invert  the  initial  file 
four  sectors  at  a time. 

The  input  cards  to  the  WAKE  program  (example  copies  are 
included  in  the  Appendix)  are  as  follows: 


Card  1:  INFLG,  N,  JOBE  (31^  format  code).  INFLG  = 0 on 

a restart,  = 1 on  a start.  N = 1 permits  a new  run  to  start 
before  completing  the  current  run;  N = 0 does  not.  JOBE  = 
value,  the  upper  minute  limit  on  current  job  execution  time. 

If  JOBE  = 0,  the  program  will  not  test  for  the  job  time. 

Card  2:  NRUNI,  CMNT  (l4,  19A4).  If  NRUNI  = value,  the 

current  run  is  given  this  run  number;  if  NRUNI  = 0,  the  run 
number  counter  in  the  common  file  is  updated  by  one.  CMNT  is 
a 19  element  vector  containing  any  desired  comments  (printed 
at  the  start  of  the  run). 

Card  3:  NSTMX,  NSTSI,  XMIN,  XMAX,  LELX,  MXHRS,  LAMIN 

(214,  3F8.3,  214).  NSTMX  sets  the  maximum  number  of  steps 
permitted  for  the  run.  NSTSI  sets  the  initial  step  value, 
typically  = 0.  XMIN  is  the  initial  x value.  XMAX  is  the 
maximum  x value  (the  program  terminates  when  reaching  XMAX). 

BE LX  is  the  initial  Ax  spacing.  MXHRS  is  the  run  time 
maximum  hours  (termination  also).  LAMIN  = 0 signals  a normal 
run  with  turbulence;  LAMIN  = 1 signals  a run  with  laminar  flow. 

Card  4:  MXRY,  LYLFF,  MXRZ,  LZLFF  (4l4).  MXRY  and  MXRZ 

are  the  maximum  number  of  mesh  points  possible  in  the  y and  z 
directions;  the  most  possible  currently  is  40.  LYLFF  and  LZLFF 
set  the  lower  boundary  flags,  = 0 implies  a free  lower  boundary 
= - 1 implies  a reflecting  lower  boundary. 

Card  5:  NRFV  (l4l4).  NRFV  (7,2)  gives  the  reflection 

properties  for  the  seven  variables  qL  , ^>,  u,  v,  w,  A,  v across 
the  two  axes.  The  first  seven  integers  give  the  variable 
properties  across  the  z axis  ( + to  - y);  the  second  seven  give 
their  properties  for  + to  - z across  the  y axis.  The  integer 
entry  is  0 if  the  variable  is  zero  at  the  axis,  and  1 if  the 
slope  of  the  variable  is  zero  across  the  axis. 
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Card  6:  VSCAV  (5F8.3).  VSCAV  is  a five  element  array 

giving  the  scaling  factors  for  q2,  p,  u,  v,  w;  typically  = 1.0. 

Card  7:  VWTFV  (5F8.3)  is  a five  element  array  giving  the 

weighting  factors  for  q2,  p,  u,  v,  w.  Typically,  weight  = 1.0 
io  given  to  q<-,  p and  u;  while  = 0.0  is  given  to  v and  w (they 
do  not  control  anything). 

Card  8:  ZEROV  ( 7 F8 . 3 ) . ZEROV  is  a seven  element  array 

giving  the  edge  values  of  q2,  p,  u,  v,  w,  A and  r, , where  q2, 
u and  u must  = 0.0,  but  nonzero  values  are  possible  for  p,  v 

and  w.  The  edge  value  for  A is  formed  internally  within  the 
program. 

Card  9:  EPSN,  EPSX,  EPSS,  ECMN,  ECMX  (5F8.3).  This  card 
gives  the  run  tolerances,  and  is  used  with  the  scaling  and 
weighting  factors  to  set  running  noise  levels  (EPSN  times  the 
appropriate  VSCAV  entry),  maximum  changes  (EPSX),  edge  toler- 
ances (EPSS),  minimum  curvature  (ECMN),  and  maximum  curvature 
(ECMX).  Typically,  EPSN  = 0.001,  EPSX  - 0.05;  EPSS  = 0.001; 
ECMN  =0.02  and  ECMX  =0.05. 

Card  10:  DXMAX,  DXMIN,  DXEMX , BUFAC,  FOUR,  DFFMN , DFFMX , 
DI'RMX  (8F8.3).  This  card  sets  the  spacing  in  the  three  direc- 
tions. in  the  marching  direction  x,  DXMAX  is  the  maximum  step 
size  permitted;  DXMIN  is  the  minimum  size,  = 0.00001;  DXFMN  is 
the  maximum  rate  at  which  Ax  can  grow  and  the  solution  march 
downstream,  = 1.5;  BUFAC  is  the  maximum  factor  (times  EPSX) 
that  a variable  may  change  before  the  step  size  is  reduced  and 
the  step  tried  again,  _ 2.0;  FOUR  multiplies  ECMN  and  ECMX  when 
curvatures  require  too  many  points,  = 1.1;  DFFMN  is  a minimum 
spacing  factor  (times  the  width  of  the  profile)  below  which  a 
point  cannot  be  inserted,  = 0.05;  DFFMX  is  a maximum  spacing 
factor  above  which  a point  must  be  added,  = 0.05;  and  DFRMX  is 
the  optimum  spacing  ratio,  traditionally  set  equal  to  2.02. 
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Card  11:  NIOLP,  NIOPP,  NFOLP,  NTOPP,  NTOLP,  NSOLP,  NIAAF, 

LOUT  (8l4).  The  first  six  parameters  control  printout  as  a 
function  of  the  number  of  steps  taken  by  the  solution.  Thus, 

If  NIOLP  = every  fifth  step  will  generate  an  intermediate 
line  printer  record  of  x,  Ax,  A , A , maximum  values  of  all 

y “ 

variables,  and  integrals  of  potential  and  kinetic  energies,  and 
momentum.  NIOPP  generates  intermediate  disk  storage;  NFOLP 
causes  a print  of  normalized  mesh  output  along  a specifically 
inputted  y and  z value;  NTOPP  generates  a disk  save  of  the 
total  variable  profiles;  NTOLP  causes  a line  print/plot  of  the 
variable  profiles;  and  NSOLP  yields  a print  of  nonnoise  turbu- 
lence to  the  line  printer.  NIAAF  sets  the  frequency  of  profile 
readjustment,  = 2 from  stability  considerations.  LOUT  = -1 
forces  only  a printer  plot;  = +1  forces  only  a numoer  plot;  = 0 
permits  both  printer  plot  and  full  profile  output  to  the  line 
printer  for  every  requested  variable. 

Card  12:  XOUTV  (10F8.3)  gives  specific  x values  at  which 

total  printouts  and  disk  storage  is  requested  (overrlces  the 
output  control  parameters  on  the  previous  card) . 

Card  13:  NSTAT,  NSTBC,  NUFF,  LIBUF,  LFBUF,  LTBUF,  NSTPR, 

IDSV  (I4l4) . NSTAT  controls  the  solution  configuration  desired: 
NSTAT  = 1 for  the  solution  of  q^  p,  u;  NSTAT  = 2 for  q^  p,  v,  w, 
v;  and  NSTAT  = 3 for  all  variables.  NSTBC  prevents  a backup  and 
retry  of  a step  for  the  initial  number  of  steps  (typically  = 8) 
set  by  NSTBC.  NUFF  sets  the  smooth  transition  to  an  XOUTV  value 
= 4.  LIBUF,  LFBUF  and  LTBUF  control  intermediate,  full,  and 
total  printouts  on  a backup  (0  or  1).  NSTPR  is  the  number  of 
steps  in  the  pressure  Iteration  loop,  < b.  IDSV  is  an  array  of 
seven  elements  giving  the  disk  storage  switch  for  each  of  the 
seven  variables  ( = 0 if  the  variable  is  not  to  be  stored,  as 
must  be  the  case  currently  for  A;  = 1 if  the  variable  is  to  be 
stored  for  future  plotting  purposes). 
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Card  14:  RE,  G,  PR,  AS,  A,  B,  BETA,  C,  S (9*8.3).  This 

card  gives  the  Reynolds  number,  (Froude  number)-",  Prandtl 
number,  and  turbulence  model  parameters  currently  set  at  AS  = 

2-5;  A = 0.75J  B = 0.125;  BETA  = 0.0;  C = 0.1;  S = 1.8. 

Card  15:  DC,  U,  RIS,  DRDZ , SCALE,  PCRIT,  SCALM,  SPORS, 

CPOR  (9F8.3).  DC  = 0.3  is  the  diffusion  coefficient;  U is  the 
free  stream  velocity  = 1.0;  RIS  = 0.25  is  the  stability  cutoff 
of  the  Richardson  number;  DRDZ  is  the  constant  background  dp0/6z. 
If  DRDZ  = 0.0,  WAKE  will  read  card  20.  SCALE  sets  the  initial 
value  of  the  scale  length.  PCRIT  = 0.005  is  the  value  of  the 
ratio  of  maximum  change  to  maximum  value  below  which  the  pressure 
iteration  is  said  to  converge  (herewith,  0.5$).  SCALM  = 0.01 
is  the  minimum  scale  length  permitted.  SPORS  is  the  square 
of  the  distance  to  the  wave-absorbing  liner,  normalized  by  A; 
and  CPOR  is  the  strength  of  the  liner  itself.  The  liner 
strength  increases  exponentially  with  a power  of  CP0R*(r2  - SPORS) 

Card  16:  SI,  S2,  S5,  S6,  S 7,  S8,  LSCAL  (8F8.3,  l4).  The 

values  of  SI  - S8  set  the  constants  in  the  dynamic  scale  equation. 
LSCAL  = 0 means  the  scale  equation  is  computed  but  does  not 
influence  the  other  variables.  Curx-ently,  s„  = S5;  s = Sb  - 
S7  = £S8.  3 4 

Card  17:  QCUT,  DIVP,  DIVF,  PNORM,  PCUT,  CMU,  XFACT,  XZERO, 

YOUT,  ZOUT  ( 10F8. 3) . QCUT  = 0.001  is  the  factor  times  the 
current  maximum  value  ol  q'  that  sets  the  computational  noise 
level  of  q . If  a local  value  of  q^  is  below  QCUT*q2  , that 
point  is  not  considered  when  constructing  the  various  integrals 
of  interest.  Re  1 is  given  the  value  of  CMU,  not  RE-1,  at 
points  beyond  a squared  distance,  normalized  by  A,  greater  than 
PCUT  (=  40) . This  procedure  is  a further  attempt  to  stabilize 
the  solution.  DIVP  = 0.1  is  the  fraction  of  DXMAX  setting  the 
minimum  value  of  the  step  size  Ax  entering  the  calculation  for 
the  corrective  divergence  effect,  the  last  term  in  eq.  (2.7). 
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When  the  program  £*new  is  below  DIVP  * DXMAX,  the  divergence 
correction  may  overwhelm  the  physics  taking  place,  typically 
during  the  first  few  steps  of  a run.  DIVF  = 2.0  is  the  accepted 
value  given  to  the  effectiveness  of  the  divergence  term,  so  that 
d/dx  ~ - u/2Axijew.  PNORM  =1.0  normalizes  the  pressure  step-size 
factors  that  yield  the  step-size  in  the  pressure  iteration  (ref.  4). 
The  pressure  step  is  At ( i ) = PNORM  * (YMAX-YMIN)  * (ZMAX-ZMIN)/PN( i) , 
where  PN  is  (10.0,  30.0,  60.0,  100.0,  150.0).  XPACT  and  XZERO 
give  the  rational  distance  downstream,  so  that  (X/D)  = XPACT  * X 
+ XZERO.  Finally,  YOUT  and  ZOUT  are  the  values  of  y and  z along 
which  the  full  printout  is  computed  for  output  display. 

Card  18:  FNAME  (8A4)  is  the  name  of  the  disk  files  to  be 

used  during  program  execution,  and  are  the  working  file  ( WAKWF) , 
the  gamma  matrix  file  (WAKGM),  the  global  intermediate  save  file 
(WAKGL)  and  the  total  profile  plot  save  file  ( WAKPL) . The  common 
file  (WAKCM)is  entered  into  WAKE  by  a data  statement. 

Card  19:  FILEN,  ISTAT,  ITYPE,  LZMAX,  LPRFL,  I FULL  ( A4,  Al, 

514).  This  card  controls  the  location  of  the  initial  profiles, 
and  provides  start  conditions  on  various  aspects  of  the  WAKE 
program.  When  ITYPE  = -1,  FILEN  is  ignored  and  the  initial 
profiles  a^e  sought  in  the  entered  profile  plot  file.  When 
ITYPE  = 0,  FILEN  gives  the  name  of  the  file  containing  the 
initial  profiles.  For  ITYFE  > 0,  these  initial  profiles  are 
assumed  to  have  been  generated  by  the  axisymmetric  program  JETMN. 
ISTAT  is  the  counterpart  of  NSTAT  and  gives  the  run  statistic 
for  the  initial  profiles  (1,  2 or  3).  When  LZMAX  = 0 the  program 
keeps  track  of  the  spread  in  Az  and  prints  and  stores  the  running 
profiles  at  the  first  point  of  nonmonotonic  behavior  in  a . 

Then  LZMAX  = 1 so  that  the  procedure  will  not  repeat.  When 
LZMAX  = - 1,  the  scales  Ay  and  Az  are  fixed  at  the  SCALM  value 
entered  on  card  15.  LPRFL  = -1  is  the  standard  operating  mode 
of  the  pressure  loop,  and  signals  to  tne  body  of  the  WAKE  program 
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that  the  iteration  is  complete.  When  LPRFL  = 0 on  initializa- 
tion, the  pressure  iterations  are  performed  prior  to  the  first 
step  I ELX  in  q^,  p,  u,  v,  w.  IFULL  = 1 forces  a total  printout 
to  the  line  printer  of  the  initial  profiles;  IPULL  = 0 does  not. 

Card  20:  ZL,  DRIZH  (2F8.3).  When  DRDZ  - 0.0  on  card  15, 

card  20  is  read.  Above  z = 0,0  and  below  z = ZL,  the  background 
density  gradient  is  se.  to  DRDZH.  For  ZL  < z < 0.0,  DRDZ  = - 10-6. 

With  the  reading  of  these  twenty  cards,  WALE  will  begin 
computation.  Beside  the  intermediate  output  described  above, 
several  integrals  and  functional  values  are  printed  every  step 
by  the  pressure  iteration  and  V/AKSC  subroutines.  The  heading 
above  each  value  is  chosen  to  give  a quick  identification  of 
that  value.  The  normalized  divergence  error  indicates  the 
closeness  to  which  our  calculation  is  maintaining  V • 7 = 0. 

The  various  kinetic  and  potential  energy  components  give  an 
indication  of  energy  distribution  within  and  outside  the  turbu- 
lent wake.  The  printout  "error"  of  the  pressure  iteration  should 
be  disregarded  since  the  pressure  iteration  scheme  currently  in 
use  has  evolved  beyond  the  usefulness  #f  "error".  Otherwise, 
printout  in  WAKOT  is  explained  when  it  is  printed.  The  plot 
file  is  stored  in  the  same  manner  that  the  initial  profile  file 
is  constructed,  only  with  mult i pie  values  of  x,  one  set  of  data 
following  the  oth^r. 

The  switches  on  the  META-4  console  permit  on-going  changes 
to  the  output  scheme  presented  above.  The  switches  used  here 
are : 

0:  terminates  a run  during  the  pressure  iteration  loop 

1:  forces  an  intermediate  line  printer  output 

2:  forces  a full  line  printer  output  at  YOUT  and  ZOUT 

3:  forces  a total  line  printer  output 

4:  forces  an  intermediate  printout  to  the  disk  file 
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5:  forces  a total  printout  to  the  disk,  file 

6:  forces  a total  line  printer  output  of  the  turbulence 

distributions 

7:  forces  a printer  picture  of  the  main  variables 

13:  terminates  a run  after  a step  in  Ax  (including  the 

pressure  iteration)  has  been  completed 
14:  prints  the  full  reflection  vector  on  run  initialization 

There  is,  of  course,  no  substitute  for  the  actual  program 
language.  A copy  of  the  FORTRAN  logic  is  presented  in  the 
Appendix.  The  potential  user  must  be  warned  that  the  weak  link 
in  the  program  is  the  capacity  of  the  pressure  iteration  loop 
to  maintain  a consistent  pressure  solution  behind  the  running 
solution  of  q , p,  u,  v and  w.  The  cross  velocities  v and  w 
will  react  rapidly  to  a degenerating  pressure  field,  but  they 
react  slowly  to  a build-up  of  the  density  field  prior  to  collapse. 

The  WAKE  program  runs  slowly  on  the  A.R.A.P,  facility;  a 24  hour 
run  is  standard  when  trying  to  extend  collapse  dynamics  to  one 
B.V.  Since  this  program  is  a "one-man"  operation,  extreme  care 
must  be  taken  in  the  input  of  the  data  cards.  Internal  data 
consistency  checks  are  simply  not  there. 
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APPENDIX 


This  Appendix  contains  the  FORTRAN  listing  of  the  complete 
WAKE  program  as  programmed  at  A.R.A.P.  The  first  listing  is 
the  mainline  WAKE  followed  by  the  listing  of  common,  followed 
by  the  called  subroutines  in  alphabetical  order.  WAKE  consumes 
all  32  K of  core  storage,  requiring  every  subroutine  called  by 
the  mainline  to  be  Localed,  and  subroutines  within  these  routines 
(particularly  WAKMY  and  WAKMZ)  to  be  Secaled.  Data  and  program 
storage  would  push  beyond  200  K on  a larger  machine. 

The  las  program  in  the  Appendix  gives  a sample  indication 
of  the  structure  of  a FORTRAN  program  generating  the  initial 
profile  file.  Following  this  program  is  a listing  of  an  example 
input  deck. 
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WAKE.S<0105) 

♦ IOCS(2501  READER « 1403  PRINTER) 

**WAKE  - STRATIFIED  SUBMARINE  WAKE  MAINLINE 
C 

c this  is  THE  stratifilo  submarine  wake  mainline 
c 

♦COPY  (CMWAK) 

c 

C INITIATE  the  run 

c 

CALL  wAKIN 
IF  (LSTFL)  20  * <:0  * 6 
C 

c integrate  the  equations 
c 

5 GO  TC  (7«bNNSs> 

6 CALL  WAKSZ 

GO  TO  ( 10  * 7 ) * NSS 

7 CALL  waksy 

GO  TO  ( 6 « 10  ) « NSS 
C 

C INITIALIZE  THE  PROFILES 

c 

e CALL  WAKPI 

C 

C COMPUTE  AUXILIARY  UUANTITIES 
C 

10  CALL  WAKAC 
C 

C CHECK  SOLUTION  VALUES 
C 

20  CALL  WAKAL 

IF  ( LMLFL ) 22*30*5 
C 

C COMPUTE  THl  SUPEREQUILIBRIUM  VALUES 
C 

22  CALL  WAKSC 

IF  (NSTAT-1)  24,30,24 
C 

C COMPUTE  PRESSURE  SOLUTION 
C 

24  CALL  WAKPC 


C 

C OU1PUT  RUN  RESULTS 
C 

30  CALL  WAKPP 

CALL  WAKOT 
IF  ( LMLFL ) 40*20,20 
C 

C AUTOPOINT  ADJUST  PROFILES 
C 

40  IF  (LIAAF)  50*20, SO 

50  ieno=nrst-*-jeno-i 

CALL  WAKAA(U2DZM,20L0V,ZNEWV,NRNZV,NRST,IEi^D*MXRZ,LZFAF,NPNZ, 
1 LZLFF  * 2 ) 

CALL  WAKAA(D2UYM,Y0LUV,YNEWV,NRNYV,NYPS,NYPE,MXRY,LYFAF ,npny, 
1 LYLFF  « 1 ) 

IF  (LFCUR)  30*60,30 
60  CALL  WAKAZ 


■ 
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A-2 


GO  TO  20 
£NL. 
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CMWmK,S(0105) 

DIMENSION  ZMV(20)  »ZV<20)  t 2PV  ( 20 ) 

C 

COMMON  IOLAY«IYPSN,NMOVE.NK. JOBE t NOUT t BUFR ( 1284 ) »BUFb(l284) 
COMMON  NWRiNWW^F,NPISIMfNCOMT,NWVEC»lBOTfITOP 
COMMON  NPOSf  JZ,  JY»LZFAF  , L YF AF , L I A AF , NMAT  ♦ N v AR « NV ART , NMK 
COMMON  LIOLF,LiOPF ,LFOLF,LTOPF,LTOLF,LSOLF,LSTFL»LMLFL 

common  ltr;mf«lburf ,lpkrf.luvcf,lovff»lfcuk, JERR 

COMMON  YM.Y, YP » FMU , XK , C V V i BBET A . HBETA » BBS . G , DX 1 1 C WS , SCAL 
COMMON  7.NEWV(40)  ,YNEWV(40)  » NRNZ  V ( 40  ) , NRNY  V ( 40  ) ,NPNZ,NPNY 
COMMON  IROWT,  IROWA,  XKOWR, 1KOWG . MOOD , IMAPV ( 4 ) « JST AT 
C 

COMMON  ZM,IYPSM,IYPEM,XMMV(6I ,XMYV(6) .XMPV(6) 

COMMON  Z,IYPS.iYPt,XZMV(b)  .XV (6) iXZPV(6) 

COMMON  ZP,IYPSP,lYPtP,XPMV<6>  »XPYV(b)  iXPPVtfe) 

COMMON  TZMV  ( 6 ) iTV<6> » TZP V ( b ) , TPYV ( 6 ) tTMYV(b) 

COMMON  A V ( fe ) , GM(lb) » KOWB ( 16 , 40 . 3 ) .K0WG(16»4Q) 

COMMON  XMAT(6)  t YMAK6)  , ZMA  I ( 6 ) , r>VEC  ( 24 ) 

COMMON  ZA,IVECA,IVECB.IRF Wfe,4> 

COMMON  SLNIL»<4)  ,GAMIU<3)  , GLO  ID  ( 3 ) , PLT1D  ( 3 ) , COMI U ( 3 ) 

C 

COMMON  NSTMX.NbTST  ,XMIN,XI“iaX,NG.NP.NSSiNRUN,KL0K(6) 

COMMON  NYPS,NYPL,MXKY,NRST» JEND»MXKZtLSCAL,LENDF»DFFMN»CMU 
COMMON  EPSN , EPbX  , EPbS  , t CMN « ECMX  » VSC  a'M  5 ) » V*  TFV  ( 5 ) 

COMMON  DXMAX,DXMIIM,DXFMX,bUFAC.DFFMX.DFRMX,DXSAV.X.XP,FCUR 

common  npts,niolp,niopp,nfolp,ntopp,ntolp,nsolp,nstpk,lprfl 

COMMON  NIAAF,NbTAl , NSTBC , MUFF , LIBUF , LFBUF , LTBUF 
COMMON  XOUTV(IO)  ,YOLUV(40)  » ZOLDV  ( 40  ) » I YPSV  ( «f  0 ) , lYpEV(40) 

COMMON  KE  i G , PR , DC  » Ai> » S , a i d f BET  A * C t U , RI S , DKjZ  i PC  R I T « SC  A«-M 
COMMON  LAMIN,MXriRS,LOUT,LZMAX,DX»CVSiD2DYM(40)  .D2qZM(40) 

COMMON  LPSXV(5) ,EPSbV(b) » tCMXV ( 5 ) i XMOM , XPE t XKE t TURBX ( 10 ) 

COMMON  FMAXV(7) tYMAXVt?) , ZMAXV ( 7 ) t TMAXV ( 7 ) ,GMAXV(7> 

CUMMOM  OEPST,UtPSI fdPORS,cPORfPNORM,PCUT,UcUT,DlVp,DIVF 
COMMON  Sl,S2,Sb,S6tS7,S8f XFACTfXZERO,SCALEtYSCAL.ZSCAL 
COMMON  ZER0V(7) ,NRFV(7,2) f IDSV ( 7 ) » LDROZ , DROZL , ORDZH 
COMMON  CMNT(19) ,FMaME<8) i L YLFF t LZLFF » YOUT » ZOUT f SP ACR ( 10 ) 

C 

EQUIVALENCE  <ZMV(1)  iZM)  , (ZV<1)  tZ)  f (ZPV(l)  »ZP> 

cart  id  0105  db  ador  4340  Db  cnt  oo4a 
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WAKAA . S ( 0 105 ) 

**WAKAA  - STRATIFIED  SUBMARINE  WAKE,  AUTOROINT  ADJUST  PROFILES 

SUBROUTINE  WAKAA (D2UM, FOLD V , FNEWV, NRNV , ISTRT ♦ IEND , MXP , LFAF*NPN« 

1 LFLFFtLL) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  DETERMINES  THE  NEW  SET  OF  POINTS 
C NEEDED  TO  SATISFY  ERROR  TOLERANCES 
C 

DIMENSION  02DM<40>  ,FOLDV<40  ) ,FNEWV  ( 40  ) , NRNv  ( 40  ) ♦ JPOS  ( 2 > 

♦COPY  (CMWAK) 

C 

OATA  JP0S/1HY«1H2/ 

C 

1000  FORMAT ( //33H  CURVATURE  TOLERANCES  RELAXED  IN  *A1«E15.5) 

C 

I OL A Y=ll 

FNEWV(l)=FOLDV(ISfRf-l) 

NRNV(1)=ISTKT-1 

DFMIN=OFFMN*(FOLDV(IENO  + l )‘-FOLOV(ISTRT-l)  ) 
OFMAX=DFFMX*(FOLDV(IEND+l)-FOLOV(ISTRT-l) ) 

istrx=istrt+lflff 

ratio=ecmx/ecmn 

FCURVrl • 0 


c 

c forward  pass  to  determine  NEW  PROFILE  POINTS 
c 

loo  JN=1 

DFNEWrl.OElO 

OU  128  J=ISTRX  , IENO 

F=FOLDV(J) 

FP=FOLDV( J+i) 

IF  (J-ISTRT)  1101*1102,1102 

1101  fm=f+f-fp 

GO  TO  llu4 

1102  FM=FOLDV( J-l) 

1104  DFM=F-FM 

DFP=FP-F 


1106 

111 

112 

113 

114 
1145 

115 

116 
117 
116 
119 
1195 


120 


DFT=DFM4DFP 
I=NRNV< JN) 

IF  (J-ISTRT)  123,1106,1106 
IF  (I)  120,120,111 
IF  (I+l-J)  120,112,112 
DFNEWrl.OElO 
IF  (JN-1)  114,114*113 
OFNEW=FNEWV< JN)-FNEWV( JN-1) 

IF  (OFT-OFMAX)  1145*1145,116 
IF  (0FT*0FT*02DM( J)*RATIO-FCURV)  115*115*116 
IF  (OFT/OFNEW-UFRrX)  126*126,116 
IF  (OFM/OFNEW-OFRMX)  117*117,119 
IF  (OFM-DFMAX)  116,118,119 
IF  (0FM*DFM*02DM(U)-FCURV)  120*120,119 
IF  (OFM-OFMIN)  120,1195,1195 
JN= JN+1 

CALL  WAKOS ( 02DM ( J-l ) , 020M ( 0 ) , OFM, DFNEW ) 

FNEWV( JN)=FM+DFM 

NRNV ( JN)=0 

JN= JN+1 

FNEWV ( JN)=F 

NRNV ( JN ) r J 
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OFNE.W  = FNEWV/<JN)-FI\EWV(  JN-l) 

IF  (QFP/DFNEW-UFRFX)  123»123«125 

123  IF  (UFP-DFMAX)  I24,124il2b 

124  IF  IUFP*0FH*U2DM(U)-FCURV)  126»126»l25 
12b  IF  (DFP-DFMIN)  126 , 1255 , 12b5 

I2b5  JiM= JN  + 1 

CALL  W AKQS ( 02DM ( J ) ,02DM(d+l) % DFP , DFNEW ) 

FNEWV ( dN)=F+OFP 
NKNV  < JN)=0 

126  IF  ldN+4-MXp)  128 « 160 i 160 
12b  continue 
JN= JN+1 

FNLWV ( JN)=FOLDV( IEND+1) 

NKNV( JN)=IEN0+1 

c bAcKWARD  pass  of  new  points  to  satisfy  Ratio  criterion 
c 

js=o 

130  JP= JN+JS 
FP=FNEWV(dN) 

FNEW V ( dP ) =FP 
NKNV (dP)=NRNV(dN) 

00  141  d=3«dN 
I=dN+2-d 
F=FNEWV( I) 

NK=NRNV(I) 

dP=dP-l 

IF  (dS)  132,132,1^1 

131  FNEWV ( dP ) =F 
NKNV(dP)=NR 

132  FM=FNEWV(I-1) 

DFP=F  P-F 
DFMsF-FM 

IF  (OFM/DFP-UFRMX)  140 , 140  » 134 
134  NKh=NRNV( 1-1) 

IF  (NR*NRM)  140,140,135 
13b  DFO=DFM 

CALL  wAK0S(D20h(NR) «D2DM(NRM) iDF0«UFP) 
fo=f-dfo 

IF  (NR-NRM-2)  137 « 136, 137 

136  NK=NR-1 
F=F0LDV(NR) 

IF  ( ABS(F-FO)-OFIVDFKMX)  138  « 138  « 137 

137  NR=0 
F=FO 

138  dP=dP-l 

IF  (dS)  140,140,139 

139  FNLWV(dP?=F 
NRNV(dP)=NR 

140  FP=F 

141  CONTINUE 
C 

C CHECK  FOR  POINTS  ttOOLD  TO  SATISFY  KATIO  CRITERION 


C 


IF 

(dS)  142,142,150 

142 

dS= 

2-dP 

IF 

( dN+dS+4-MXP ) 143«143»160 

143 

IF 

(dS)  150,150,130 

C 

l. 
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c check  whether  profile,  must  be  reorganized 
c 

150  LFAF=0 
NR=ISTRT-1 
JN*JN+JS 

00  156  J=1»JN 
NKNrNRNV ( J ) 

IF  (NRN-NR)  151 *155*151 

151  LFAF=1 

GO  TO  156b 

155  NR=NR+1 
15b  CONTINUE 
15&5  NHN=JN 

IF  (FCURV-1.0)  15«»158,157 
157  WRITE  ( NOUT  « 1000 ) JPOS(LL) tFCURV 

156  RETURN 
C 

C MAXIMUM  POINTS  EXCEEOEU  ON  AUTOPOINT  ADJUSTMENT 
C 

160  FCURV=FCURV*FCUR 

IF  (FcURV-l.OElO)  lOO.iaOflSO 
180  LFCUR=LL 

return 

ENU 
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♦*WAKaC  - stkatifieo  submarine 

SUBROUTINE  UAKAC 


make « auxiliary  calculations 


this  SUBROUTINE  IN  the  wake  program  performs  a number  of 

AUXILIARY  TO  The  SOLUTION  OF  THE  SYSTEM  OF  EQUATIONS 

1)  COMPUTES  MAXIMUM  ABSOLUTE  VALUE  FOR  EaCH  VARIABLE  AND 
LOCATIONS  OF  SmME • MAXIMUM  STEP  CHANGE  OF  VARIABLES, 
MAXIMUM  X DERIVATIVES 


Tasks 


y and 
and 


COMPUTES  FOR  EmCH 
WITH  RESPECT  Tu  2 

COMPUTES  FOR  EaCH 
WITH  RESPECT  TU  Y 


2 ROW  THE  MAXIMUM  SECOND  DERIVATIVE 
FOR  ALL  Y 

Y COLUMN  THE  MAXIMUM  SECOND  DERIVATIVE 
FUR  ALL  2 


♦ COPY 


CONSTRUCTS  THE  MESH  VALUE  PRINTOUT 
COMPUTES  THE  APPROPRIATE  LENGTH  SCALES 

dimension  value (40) , total < 7,40,2*  ,psivuo) 

(CMwAK) 

EQUIVALENCE  ( Tut aL (1,1,1) ,K0WG(ltl> ) , (VaLUE(I) ,NRNZV(2) ) 
EQUIVALENCE  (PSIV(1 > «2NEWV(1) ) , ( GM ( 1 ) , TKE ) , (GM(2) ,VKL> 
EQUIVALENCE  (GM(3) ,KPE) , (GM(4) ,wPE) , (GM(5) ,WKE) 
EQUIVALENCE  (GM(b) ,AREA) , (GM(7)  *PSIM) , (GM(8) ,XLIFT) 


2EKO  PERTINENT  VARIABLES 

IULAY=5 

MOOD=1 

lfl=nptsn-nstst 

PLANE=1»0/FLOAI (LYLFF+2)/FL0aT (L2LFF+2) 

XMuM=0.0 

XPE=0.0 

TKE=0.0 

VKE=0.0 

RPE=0.0 

WHL=0.U 

WKE=0.0 

AREA=0.0 

PSIM=0.0 

XLIFTsO.O 

CALL  SFVFL(0.0,PSIV,MXRY) 

CALL  SFVFL(0.0,FMAXV,4^NVART) 

DO  40  1=1,2 
DO  35  UY=1,4Q 

CALL  SFVMV(ZERUV,101AL(1,UY,I),NVART) 

CONTINUE 

CONTINUE 

CALL  SFVFL ( 0 • 0 , D2UYM , 40 ) 

CALL  SFVfL(0.0,D2U2M,40 ) 

CALL  WAKSE(ECMX,ECMXV) 

nypsx=nyps+lylff 
00  50  uy=nypsx,nype 

IF  ( YOUT-YOLUV( JY) ) 60,55,45 
IF  ( YOUT-YOLDV( jY+1 ) ) 55,50,50 
CONTINUE 
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L 

I 


I 


55  RATY=< YOuT-YOLOV(UY) )/( YOLOV( JY+1)-Y0LDV(UY) ) 

C 

C INITIALIZE  FOR  PASSING  THROUGH  PROFILE 
C 

60  IR0WAz2 

IR0WR=4 

IF  (LSTFL)  96 ,96,95 

95  IR0WRz2 

96  NRSTXrNRST+LZLFF 
NRENDzNRST+jEND-1 
nhos=nkst-i 

CALL  WAKRR(NKST-1«ZV) 

CALL  WAKRR(NRST,ZPV) 

DO  250  NRsNRSTX.NHEND 
PSiHrO.O 
C 

C REaD  three  SURROUNOING  ROWS  AND  test  for  selected  Z in  DOMAIN 

c 

IF  (NR-NRST)  104,103,103 
103  NPOS=MR 

CALL  SPVMV(Z«ZM,NWWZF ) 

CALL  WAKMR (2,1) 

CALL  SFVMV(ZPtZ,NWWZF) 

CALL  WAKMR ( 3 « 2 ) 

CALL  WAKRR(NR+1»ZPV) 

GO  TO  105 
10**  ZM=Z+Z-ZP 
105  DZM=Z-ZM 
DZPzZP-Z 
DZT=ZP-ZM 
LYFAFrO 

IF  (ZOUT-ZM)  107  « 1065 , 106 
10t>  IF  (ZOUT-Z)  1065.107,107 
1065  LYFAFsl 

RATZ=(ZOUT-ZM)/iZ-ZM) 

C 

c ste.p  through  all  y points  computing  auxiliary  quantities 
c 

107  YrYOLDV(lYPS-l) 

YP=YOLOV( IYPS) 

CALL  WAKMP(IROWR,NR-*l,  I YPS-1 , XP YV , 1 ) 

CALL  WAKMP(IROwR,NR,IYPS-1,XV,1) 

IF  (NR-NRST)  100,109,109 
10b  CALL  WAKRF ( XPYV , XMYV , 2 ) 

GO  TO  110 

109  CALL  WAKMP(IR0WR,NR-1,IYPS-1,XMYV,1) 

110  IYPSXzIYPS-fLYLFF 

IF  (LFL)  112,112,111 

111  CALL  WAKMP(IR0WA,NR,IYPS-1,AV,1) 

112  CALL  WAKMP(IR0WR,NR+1,IYPS,XPPV»1) 

CALL  WAKMP(IR0WR,NR,IYPS,XZPV,1) 

IF  (NR-NRST)  113,119,114 

113  CALL  WAKRF(XPPV,XMPV,2) 

GO  TO  115 

114  CALL  WAKMP(IROwR,NR-1, IYPSiXMPV.l) 

115  DO  200  JY=IYPSX,IYPE 

IF  (JY-IYPS)  116,117,117 

116  CALL  WAKRF(XPPV,XPMV, 1) 

call  wakrf(xzpv,xzmv,d 
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U 

IJ 
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CALL  WAKRF(XMPv,XMMV,l) 

ym=y+y-tp 

GO  TO  130 

H7  IF  (LFL)  120,120,119 
H9  CALL  WAKMP(IR0wA,NR«UY«AV,1) 

120  CALL  SFVMV(XPYV,XPMV,NMOVE) 

CALL  SF VMV ( X V » XZMV , NMQVE ) 

CALL  SF  VMV  ( XMY  V » XI*MV  » NMOVE  ) 

YM=Y 
Y = YP 

YP=Y0LUV( JY+1) 

CALL  WAKMP(  IR0*R,NR+1,JY+1.XPPV,1> 

CALL  WAKMP<  iROwR.NR.  JY+1 . X*PV , 1 ) 

IF  (NR-NRST)  125 « 120  <126 
12b  CALL  WAKRF(XPPV,XrPV«2) 

GO  TO  130 

12b  CALL  wAKMP(IROwR«NR-1»JY+1«XMPV.1> 

130  DYM=Y-YM 

dyp=yp-y 

dyt=yp-ym 

c 

C PthFORM  AUXILIARY  COMPUTATIONS  FOR  EACH  VARIABLE 
C 

DO  145  I = 1 « NwVLC 
IF  (LFL)  1406, 140b, 1402 

1402  TEM=ABS(XV(I)-AV(1) ) 

IF  (TEM-TMAXV(l))  140b, 140b, 1403 

1403  TMAXV ( I ) =TEM 
140b  TEM=ABS(XV(I) ) 

IF  (TEM-FMAXV(l))  142 « 142 « 1407 
1407  FMAXV  ( I ) =TEM 
YMAXV ( I ) =Y 
ZMAXV ( I ) =2 

IF  (1-1)  1406.1408,142 


1408 

142 

143 

144 


144b 

1447 

1446 

14b 


nkl=nr 


JYL=JY 

IF  (LIAAF)  I43.l4b.143 
IF  (I-NVAR)  144,144,145 

D2DY=( (XZPV(I)-XV(D )/DYP-(XV(I)-XZMV(I) )/OYM)/QYT/ELMXV( 

D2DZ= ( (XRYV(I)-XV(I) )/OZP-(XV(I ) -XMYV ( I ) )/QZM) /d2T/ECMXV( 


I) 

I) 


TEM=ABS(D2DY) 

IF  (TE«-02DYM(UY) ) 1447,1447,1446 

D2DYM(JY)=TEM 

TEM=ABS(D20Z> 

IF  (TEM-D2DZM(NR) ) 145, 14b, 1446 

D2DZM(NK)=TEM 

CONTINUE 


C 

C ST0Rt  MESH  OUTPUT  VALUtS 
C 

IF  (LYFAF)  lblo, lblb, 150 

150  DO  151  1=1 , NWVEC  w 

TOTAL ( I , JY,1)=XMYV ( I ) +RATZ* ( XV ( I ) -XMYV ( I ) ) 

151  CONTINUE 

1515  IF  (YOUT-YM)  155,153,152 

15*  IF  (YOUT-Y)  153,135,155 

153  DO  154  1=1, NWVEC  „ 

TOTAL(I,NK,2)=XZMV(I)+KATY*(XV(I)-XZMV(in 

154  CONTII^UE 


- 


I 
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155  IF  (JY-JYL)  158,156*156 

156  VALUE ( NR ) zXV ( 1 ) 

C 

c computation  OF  integhals 
C 

158  CALL  WAKLL(TEM,l) 

IF  (TEM)  200*159*200 

159  SUMF=PlANE*DYT*OZI 
LZFAFrQ 

IF  (NR-NRST)  1591*1592,1592 

1591  SUMF=0.5*SUMF 
LZFAFrl 

1592  IF  (JY-NYPS)  1593* 1594, 1594 

1593  bUMF=0.5*SUNF 
LZFaF=1 

1594  XM0M=XM0M+XV(3)*(U+XV(3) )*SUMF 
TKE=TKt+XV(l)*(U+XV(3) )*SUMF 
TEM=(XV(4)-ZER0V(H) )**2+CXV(5)-ZER0V(5) )**2 
TEm=(xV(3)*XV(3)+1EM)*(U+XV(3)  )*SUMF 
RPL=RPL+X V ( 2 ) *XV ( 2 ) *SUMF 

VKE=VKL+TEM 

XLIFT=XLIFT+(XV(5)-ZEROV(5) )*SUMF 
IF  (XV(l)-acUT*FMAXV(l) ) 162,162,161 
161  XKt=Xpt+XV(2)*Z*SUMF 
area=akea+sumf 
GO  TO  165 
lfa2  WKL=WKt+TEM 

TLM=XV(2)*XV(2)*SUMF 

XPE=XPt+TEM/2.0 

WPtrWPE+TEM 

165  IF  (LZFAF)  200,166,200 

166  PSIH=PSIH+0.5*OYM*(XV(5)+X*MV(5)-2.0*ZEROV(S) ) 

PS1V(  JY)=PS1V(  JY)+0.5*DZM*(X\/(4)  + XMYV(4)-2,0*ZEROv(4)  ) 
TEM=0,5*(PSlV( JY)-P51H) 

IF  (ABS(TEM)-ABS(PSIM) ) 200,200,167 

167  PSlMrTEM 
200  CONTINUE 
250  CONTINUE 

TKE=TK£/2.0 
VKE=VKE/2.0 
WKL=WKL/2,0 
XKE=TKt+VKE 
XLlFTr-XLIFT 
RPL=G*RPE/2.0 
WPL=G*WPE/2, 0 
XPE=G*XPE 
C 

C COMPUTATION  OF  GROSS  SCALE  LaMBOAS 
C 

NPOS=NKL 

CALL  WAKRR(NRL,ZV) 

YrYULOV ( JYL) 

XV ( 1 ) rFMAXV ( 1 ) 

TtM=-XV(l  )/4,0 

IYPSX=JYL+1 

LFL=0 

UO  260  JY=IYPSX,IYPL 
XMY V ( 1 ) =XV ( 1 ) 

YM=Y 


m 


m 


mi 


mm 


mm 


mm 
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CALL  UAKMP(IKOwR'NRL'JY'XV'l) 

y=yoldv< jy) 

IF  (TEM-XVU))  256,2t>2.252 
If  (LFU  260 '254*260 

YP=YM+(  Y-YM)*(XMYVa)-TE«)/(XMYV(l)-XV(l)  ) 

LFL=1 
60  TO  260 
LFL=0 
CONTINUE 

EPSSV(3)=2.0*C*(YP-Y0UT ) 

NKSTX=IV‘>L+1 

LFL  = 0 

00  260  NR=NRSTX,I\.K£NU 
IF  (TEM-VALUE(NR) ) 276*272*272 
IF  (LFL)  2 S 0 *274  *280 

ZP=Z0L0V  (NR-1 ) ♦ ,(2‘0LUV  (NR  ) -ZOLDV  ( NR-1 ) )* (VALUE (NR* 1 ) - TEN ) / 
1 ( VALUE(NK-l J-VALUE(NR) ) 

LFL=1 
GO  TO  260 
LFL  = 0 
CONI INUE 

EPSSV(1)=2.G*C*(ZP-Z0UT) 

return 

ENU 

ID  0105  0B  A00R  4390  DB  CNT  0206 
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WAKaL.S(0105) 

**WAKaL  - STKATIFIEO  SUBMARINE  WAKE.  AUXILIARY  LOCP  FOR  EXECUTION 

subroutine  wakal 

c 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  EXAMINES  SOLUTION  VALUES 
C 

♦COPY  (CMWAK) 

C 

DATA  JMSGF, JMS&R.UMSGB. JMSGS/2HST.2HRS.2HBU.2H0K/ 

C 

1000  FORMAT (//47H  FULL  OUTPUT  AT  FIRST  LOCAL  MAXIMUM  OF  2 LAMBOA) 

1001  FORMAT ( //35H  UX  LtSS  THaN  MINIMUM  ALLOWEQ  VALUE) 

C 

c check  run  position 
c 

10LAY  = 6 

IK  (LMLFL)  100, 100,122 

100  XP=X 
NPTSN=NPTS 

IF  (LSIFL)  102,101,102 

101  NPTSN=NPTS+l 
XP=X+DX 

102  CALL  WAKFS(NIAAF,LIAAF) 

LI OLF  = 0 

CALL  WAKFS(NI0LP,LIULF) 

CALL  WAKFS(NIOPP.LIUPF) 

CALL  WAKFS(NFOLP,LFOLF) 

CALL  WAKFS(NTOPP,LTUHF) 

TALL  WAKFSlNlULP.LTULF) 

CALL  WAKFS(NSOLP,LSULF) 

JSTATsJMSGF 

IF  ( LSTFL)  103,110,105 

103  IF  (LPRFL)  104,170,170 

104  LiOLFrl 
LIOPF  = 0 
LFOLFsO 
LTULF=0 

ltopf=o 

LSULFsO 
JSTAT=JMSGR 
GO  TO  170 
10b  LIOLF=l 
LFOLF=l 
LTOLFrI YPSN 
LSOLF=0 
LIOPF=NIOPP 
LTOPFrNTOPP 
LIAAFsNl AAF 
GU  TO  1&4 
C 

C SET  UP  FOR  NEXT  STEP 
C 

HO  CALL  WAKSE ( EPSS, EPSSV ) 

CALL  WAKSE ( EPSX , EPSXV ) 

JSTAT= JMSGS 
DXSAVrDX 

DO  11H  IXBRK=1,10 
XdRK=XOUTV(IXBRK) 

IF  (X-X8RK)  112,114,114 


mumhm 





DO  1122  1 = 1 , NUFF 
TEM=FLOAY(I) 

IF  ( X+TEM*DX*XbRK ) 1122,1124,1124 
2 CONTINUE 
GO  TO  114 

4 DX=(X0RK-X)/TEM 

IF  (1-1)  1126 , 113 » 1126 
b XP=X+OX 
GO  TO  116 
i XPsXBRK 
ib  LIOLFsl 
LFOLFsl 
LlOLFrl 
LSULFsl 
LlGPFsNIOPP 

ltcpf =nt  opp 

IF  (LTRNF)  119* 116 , 119 
i CONTINUE 

» IF  (XP-XMAX)  116,117,117 

' LTRNFrl 
LENOFsl 
GO  TO  1135 

) IF  (NpTSN-NSTMX)  119,117*117 

) UXI=2.0/OX 

INITIALIZE  for  current  integration  STEP 


LMLFLsl 
IF  (LBURF) 
NSS=3-NSS 
RETURN 
LMLFLsO 
IF  (LPKRF) 


120,12U,121 


190  * 124 , 190 


error  check 


1241 


1243 


1244 

C 


CALL  DVCHK(lSW) 

GO  TO  (1241,1242) *ISW 

LUVCF=1 

LIOLFsl 

CALL  O'/ERF  ( ISW  ) 

GO  TO  (1243,1244,1244) , ISW 

LOVFF=l 

LIOLFsl 

IF  (LOVCF+LOVFF)  17b, 130*175 


check  for  backup  anu  determine  new  dx 


130  LBURFsO 

DO  136  X = 1 * NVAK 

IF  (TMAXV(I)-EPSXV(I)*BUFAC)  136,136,135 
13b  IF  (NPTSN-NSTBC ) 136,136,1355 
1355  LBURFsI 
136  CONTINUE 

OXNsOXSAV*  AflS  ( OXF1"'  X ) 

DO  162  Isl , NVAR 
VMAXsGMAXV(I) 

IF  (LBURF)  144,142,144 
142  TEMsFMAXV(I) 


A-l4 


IF  (TEM-VMAXI  1435,1435.143 

143  GMAXV(I)=T£M 
143d  VMAX=TEM 

144  IF  (VWTFV(I))  162*162,145 

145  IF  ( VMAX-VSCAV( I )*EFSN)  1455*146*146 
1455  VNAX=VSCAV( I )*EPSN 

146  AERK=ABS(EPSX) 

IF  ( EPSX ) 147 « 146 « 148 

147  AERR=AEKR*VSCAV(I> 

GO  TO  149 

14b  AERK=AEKR*VHAX 
149  AERK=ALKR/VwTFV(I  ) 

TlM=TKAXv(I J/OX 
IF  (TEM)  151*162*151 
151  OXX=AERR/T£M 

IF  (OXX)  162*162,160 

160  IF  (QXX-DXN)  161,162,162 

161  OXN=OXX 

162  CONTINUE 
DX=OXN 

IF  (DX-UXI1IN)  1620,1621,1621 

1620  LIOLFzl 
LTRNFsl 

WRITE  ( NOUT , 1 OU 1 ) 

GO  Tu  190 

1621  IF  (DXFMX)  1&3 , 1623 « 1622 

1622  IF  (OX-DXMAX)  163*163,1623 

1623  DX=DXNlAX 

163  IF  (LBUKF)  163U . 1&4 « 1630 

1630  JSTATrdMSGB 
LTRNFrl 

IF  (UXSAV-DX*DXFMX)  1631*1632*1632 

1631  OX=OXSAV/uXFMX 

1632  LIOLFrO 
LFOLF=0 
LT0LF=0 
LSQlFzO 
LlOPFrO 
LTOPFrO 

IF  (LIBUF)  1633,1634,1633 

1633  L10LF=1 

1634  IF  (LFBUF)  l63D,l63b«1635 

1635  LFOLFrl 

1636  IF  (LTBUF)  1637,1638*1637 

1637  LTOLF=l 

1636  IF  (LIBUF+LFBUF+LIBUF ) 190*101*190 
C 

C SUCCESSFUL  INTEGRATION  STEP 
C 

164  IF  (LZMAX)  1646,1642*165 

1642  IF  (EPSSV(1 J-ZSCAL)  1644*165.165 
1644  L2MAX=1 
LIOLFsl 
LFOLFsl 
LTOLFrl 
LI0PF=NI0PP 
LTOPFrNTOPP 
WRITE  ( NOUT , 1000  ) 

GO  TO  165 
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1646 

16b 

170 

17b 


19  U 
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EPSSV(1)=SCALM 
LPSSV ( i ) =SC ALM 
ZSCAL=EPSSV<1> 

YSCAL=EPSSV<3) 

CALL  SFV(*1V(FMAXV.GMAXV,NWVLC) 
LMLFL=-1 

CALL  wAKCS(l«LiOLF ) 

CALL  WAKCS  ( '£ « LFOLF 
CALL  WAKcS ( 3 « LI OLF ) 

CALL  WAKCS ( 4 « L10PF ) 

CALL  WAKCS ( b « LTOPF ) 

CALL  WAKCS(6«LS0LF ) 

CALL  wAKCS(U,LTRNF ) 

CALL  WAKCS<l3.LTRNF) 

RETURN 

LNU 
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WAKAY , S ( 0 105 ) 

♦♦WAKAY  - STRATIFIED  SUBMARINE  WAKE,  autopoint  adjust  y 
SUBROUTINE  WAKAY(NRX) 

c 

C THIS  SUoROUTINE  IN  THE  WAKE  PROGRAM  REORGANISES  A ROW 
C 

DIMENSION  IYPNVUO) 

♦COPY  (CMWAK) 

EQUIVALENCE  ( YOLOV ( 1 ) , 1YPNV < 2 ) ) 

C 

C MAKE  WORKING  COPY  OF  NRNYV 
C 

NPOS=nRX 

CALL  WAKRH(NRX,2V) 

IF  (IYP5»)  160*160*100 

100  IYS=0 

DO  108  I Y = 1 » NPNY 
I=NRNY V ( I Y ) 

IF  (I)  10b, 106,101 

101  IF  (IYS)  102,102*104 

10 2 IF  ( I-IYPS ) 104,103*103 

103  ITS=1Y 

104  IF  (IYPE-I)  10b, 105*105 

105  ite=iy 

10b  I YPNV ( I Y ) = I 

106  CONTINUE 
C 

C CHECK  FOR  MIO-POINT  ADUED  ON  EITHER  END 
C 


111 

IF 

IF 

(IYS-1)  114*114,111 

(NRNYV( IYSJ-IYFS)  114*112*114 

112 

IF 

<NRNYv< IYS-1)  ) 113,113,114 

115 

ITS 

=IYS-1 

114 

IF 

(IYE-NPNY)  115*120,120 

115 

IF 

(NRNYV(IYE)-IYPEI  120*116*120 

lib 

IF 

(NPNYV( IYE+1) ) 117,117,120 

117 

I YE 

=1 Yt+l 

C 

C CHECK  WHETHER  NRX  ROW  IS  VOIDED  BY  ADJUSTMENI 

c 

120  IF  (IYS-IYE)  121,121*1201 
1201  IF  (IYE-I)  121,121,1203 
1203  IYS=IYE 

C 

c forward  reorganization  pass 
c 

121  JYMsO 
IYS=MMAX(1,IYS-1) 

I YE=MMIN(NPnY» IYE  + 1 ) 

00  150  IY=IYS,1YE 
JTT=NYPS+IY-2 
JYF=IYPNV( IY) 

1 = 1 Y 

122  1=1+1 

IF  (I-NPNY)  123,123*124 
125  JYP=IYPNV(I) 

IF  (JyP)  122*122,125 

124  JYP=loOOO 

125  IF  (JYF)  126*126*132 
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IF  (JYT-JYM)  150, 150*127 
IF  (JYT-JYP)  12b«150*150 
CALL  WAKIY( JYM* JY1 ,JYP*1Y) 

GO  TO  148 

IF  (JYT-JYF)  133  # 148  * 142 
, IF  (JYT-JYM)  134*134* 145 

jym=jyf 
GO  TO  150 

1 IF  (JYT-JYP)  145*150*150 

> CALL  SfVMV (ROWB( 1 « JYF *2) *ROWB ( 1* JYT»2 ) *NM0VE) 
) IYPI\IV(IY)=-JYT 
JYM=JYT 
) CONTINUE 

backward  reorganisation  pass 


DO  154  JY=1YS«1YE 
jYT=NYPS+lYS+lY£.-JY-2 
lY=lYS+lYt-JY 
JYF  = IYPI\IV(IY) 

IF  (JYF)  154*151,152 

151  JYM=IABS(IYPNV(IY-1) ) 

JYP=IABS( IYPNV(IY+1) ) 

CALL  WAKIY( JYM,JY1 ,JTP, tY) 

GO  TO  153 

152  CALL  SFVMV (ROWb ( 1 « JYF ,2 ) ,ROWB ( 1 ♦ JYT ,2 ) ,NMOVE ) 

153  I YP1MV  ( I Y ) =- JYT 

154  CONTINUE 
IYP5=NYPS+IYS-1 
IYPE=NYPS+IYE-3 

160  CALL  WAKvJR(NRX,zV  ) 

RETURN 

END 

CaRT  10  0105  OB  AUDR  32D0  OB  CNT  00B6 
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WAKAZ , S ( 0105 ) 

**UaKaZ  - STRATIFIED  SUBMARINE  WAKE,  AUTOPOINT  ADJUST  Z 
SUBROUTINE  WAKAZ 
C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  REORGANIZES  THE  SOLUTION  FILE 
C 

♦COPY  (CMwAK) 

C 

C check  profile  SHIFT 
C 

I0lAY=12 

MOOD=-l 

IF  (L2LFF)  101*102,102 

101  NKSTN=2 
GO  TO  106 

102  NRSTN=(MXKZ-NPN2)/2+2 

IF  (NRST  + ( JENO-1  )/2-MXR2/2)  107*106,107 
107  LZFAFzl 

106  IF  (LYLFF)  111*112,112 

111  NYPSN=2 
GO  TO  116 

112  NYPSN= (MxRY-NPNY ) /2+2 

IF  (NYHS>+(NyPE-NYPS)/2-MXRY/2)  117*116,117 
117  LYFAFrl 

H6  IF  (LYFAF)  1181,119*1181 
C 

C SOLUTION  PROFILE  MUST  BE  REORGANIZED  IN  Y DIRECTION 
C 

H8l  NYPS=NYPSN 

NYPE=NYPS+NPNY-3 
DO  1185  JZ=1«NPNZ 
NHF=NRNZV< JZ) 

IF  (NRF)  1165,1165,1182 
1162  CALL  WAKaY(NRF) 

1166  CONTINUE 

UO  1168  1=1 , NPNY 

JZ=NYPS+I-2 

YOLUV ( JZ ) = YNtWv ( I ) 

1188  CONTINUE 

119  IF  (LZFAF)  120,160,120 
C 

C SOLUTION  PROFILE  MuSl  BE  REORGANIZED  IN  Z DIRECTION 
C 

120  NRST=NRSTN 
JcNU=NPNZ-2 

C 

c forward  reorganization  pass 
c 

NRK=0 

DO  150  JZ=1,NPNZ 
NRT=NRSWJZ-2 
NKpsNRNZv ( JZ ) 

J=JZ 

122  J=J+1 

IF  (J-NPNZ)  125,123,124 
125  NRpsNRNZ V ( J ) 

IF  (NRP)  122*122,125 
12**  NRP=10000 

125  IF  (NRF)  12b, 126, 132 


I 

Q 

n 

u 

u 

y 

LI 

y 

y 

y 

Ll 

o 

y 

D 

y 

o 


12b 

127 

126 

132 

133 
13*+ 

142 

14b 


146 


IF  (NRT-NKM)  150*150*  127 
IK  (NRT-NKP)  126*150* 150 
CALL  WAKIZ(NRM.NK1 ,NKP, JZ) 
GO  TO  148 


IF  (NRT-NRF) 

IF  (NRT-NKM) 

nrm=nrf 
GO  70  150 
j;F  (NRT-NKP) 
fJP05=NRF 

CALL  WAKKK(I\|KF«ZV) 
npos=nRT 

CALL  WAKwR(NRT,ZV) 
NKN^V( JZ)=-NKT 
NRM=NRT 
CONTINUE 


133*146*142 

154*134*145 


145*150*150 


150 
C 

C tsALKwARD  REORGANIZATION  PaSS 
C 

UU  154  UZ=1«NPNZ 
NKTrNRST+NPNZ-jZ-1 
IZ=NPN^+1-JZ 
NKF=NRN^V( IZ) 

IF  (NRF)  154*151,152 

151  NKM=IABS(NRN2V(1Z-1) ) 
NKP=IABS(NRN^V(IZ+1) ) 

CALL  WAKIZ(NKW,NKT ,NKP, IZ) 

GO  TO  153 

152  NPOSrNKF 

CALL  LAKRR(NKF,ZV> 

NPOSrNKT 

C Al.L  WAKWR(NRT,ZV) 

153  NKNZV(I2)=-NRT 

154  CONTINUE 

DO  156  I=1«NPNZ 

JZ=NRsT+I-2 

ZULDV(JZ)=ZNEWV(I) 

15b  CONTINUE 
160  RETURN 
END 

cart  id  0105  db  audr  21-70  db  cnt  00C4 
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WAKCL.S(0105) 

**WAKcL  - stratified  submarine  wake,  compute  length  scales 

SUBROUTINE  WAKCL(XXVfFTEM) 

C 

C THIS  SUBROUTINE  IN  THE  WAKL  PROGRAM  COMPUTES  THE  VERTICAL  AND 
C SCftLE  LENGTHS  bT  COMPARISON  WITH  THE  RICHARDSON  LENGTH 
C 

DIMENSION  XXV(b) 

♦COPY  (CMWAK) 

c 

SCALrZSCAL 

IF  (FTEM)  103»140«l03 
103  TLM=SQKT(RIS^ABS(XXV(1)/FTlM)/G) 

IF  (SCAL-TEM)  140  » 14U«11Q 
110  scal=tem 

IF  (SCAL.SCALM)  120*130,130 
120  SCAL=SCALM 

130  SCAL=SQRT< AbS(SCAL^^SCAL) ) 

140  IF  (LsLAL)  130,160,150 
150  SCALE=XX V ( 6 ) 

RETURN 

160  SCALE=2«0*YSLAL#YSCAL*SCAL/( YSCAL^YSCAL+SCaL^SCaL) 
return 
END 

cart  id  0105  db  audk  3A90  ub  cnt  002E 
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WAKCS.S(0105> 

**WAKCS  - STRATIFIED  SUBMARINE  WAKE,  CHECK  DATA  SWITCH 
SUBROUTINE.  WAKCS  ( LL  ♦ LFL ) 


C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  CHECKS  THE  LL 
C FOR  POSSIBLE  OUTPUT  TO  LINE  PRINTER  OR  DISK  FILE 
C 

♦COPY  (CMWAK) 

C 


DATA  SWITCH 


CALL  OATSW ( LL « 1 SW ) 

GO  TO  (10*20) «iSW 
10  LFL=1 

LlOLFrl 
20  RETURN 

END 

CaRT  10  0105  OB  aDDR  2OE0  db  LNT  OOiC 
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WA*UG.S(0105> 

**HAKdG  - STHATIFIEO  SUBMARINE  WAKE,  DETERMINE  DENSITY  GRADIENT 
SUBROUTINE  WAKUGUPUS) 

C 

C THIS  SUBROUTINE  IN  TH£  WaKE  PROGRAM  DETERMINES*  THE  BACKGROUND 
C DENSITY  GRADIENT  WHERE  APPLICABLE 
C 

♦COPY  (CMWAK) 

C 

IF  (LOROZ)  I0,i00»l0 
10  IF  (ZPOS)  20*30,30 

20  IF  (ZP0S-DKD2L)  30,30,25 

2b  DRDZ=-0, 000001 

RETURN 

3o  DKU2=DRDZH 

100  RETURN 
ENu 

CART  iU  010b  Dd  AUDK  i+bAO  DB  CNT  0020 
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WAKQS.S<0105) 

♦♦WAKuS  - STRATIFIED  submarine  wake,  obtain  delta  spacing  ratio 

SUBROUTINE  WttKDS(CUKVM,CUKVP,DF,DFO) 


THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  COMPUTES  a NEW  OF  AS  SOME 
FRACTION  of  the  Old  SO  AS  TO  oALANCE  the  three-point  profile 


♦COPY  (CMWAK) 

c 

DFRM=0»99*DFAMX 
IF  (CURVP)  190,191,190 

190  l£MB=SUKT<CURVtt/CURVP) 

IF  (TEMB-OFRM)  192,192,191 

191  TEMB=OFRM 
GO  TO  194 

192  TEMA=l • 0/DFRM 

IF  (TEMB-TEMA)  193,194,194 

193  TEMd=TEMA 

194  DF=uF/(  1,0+TEMtj ) 

IF  (OF/OFO-DFKM)  196,196,195 

195  DF=0F0*UFRM 

196  RtTURN 
ENU 

CART  ID  0105  OB  AODR  2?00  OB  CNT  002C 


A-24 

WaKEC.S(0105) 

**WAKEC  - STRATIFIED  SUBMARINE  WAKE,  EDGE  CONDITION  CHECK 
SUBROUTINE  WAKEC(LFL) 

r 


C THIS  SUBROUTINE  IN  THE  WAKE 

c and  checks  for  satisfaction 
c 

♦copy  ( cmmah ) 


PROGRAM  turns  the  IMPLICIT  UPSWEEP 
OF  THE  BOUNDARY  EDGE  CONDITION 


c 


lfl=i 

CALL  SFVMV<  aVNXMaI  ,NWVLC) 

DU  100  I = 1 « N W V £ L 
XMAT(I»  = XMAT(I)-GMU)*^£ROW(l> 
loo  CONTINUE 

DU  HO  1:1 , NVAR 

IF  (ABS(XMAT(I)-2EROV(I))-£PSSV(I))  110,110,120 
110  CONTINUE 

call  sfvmvmxmat, av,nwvec) 

CALL  WAKMP( IKOwa,NR«UY, aV,2) 

return 

120  LFL=0 


return 

ENU 

CaHT  ID  0105  DB  aDOR  4JiO  06  CNT  002C 
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WAftFS .S ( 0105 ) 

**WAKpS  - STKATIPIE.0  SUBMARINE  WAKE,  FLAG  SET 
SUBROUTINE  WAKFS(NPP.LFL) 

C 

C THIS  bUBROUTlNE  Iw  THE  WAKE  PROGRAM  SETS  THE  OUTPUT  FLAGS 
C 

♦COPY  (CMwaK) 

L 

LFL=0 

IF  (NPP)  10k, 102, 100 

100  IF  (MOU(NPTSN,NPP) ) 102,101,102 

101  LFL=1 
LIOLFsl 

10k  KLTUKfM 
END 

CaKT  ID  0i05  OB  AOOR  2^00  Dd  CNT  OOlC 
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WAK1N.S(0105) 

**k4KI^=F;^rM4HINE  "4KE' INPUT  C4RDS  AN0  STal,T  RUN 

c 

c this  subroutine  in  the  wake  program  starts  the  run 
c 

DIMENSION  ETIDls) « LFLAG ( b ) , COMN(2) *FILEN(2) .FILID(3) 
•CUPY  IOflT-“*J-J^HV<5».FT1D(3..HSCALfJ:;FILID‘3> 

EQUIVALENCE  ( LTRNF , LFLAG (1)),(AV(1),FILID<1)) 


c 


data 

Data 

data 


c 

1000 

1001 
1002 

1003 

1004 

1005 
I02o 
2000 
2002 
2006 


C0MN/4HWAKCi1HM/, JERKX/2H  / i BLANK/4H 
HSCAL/3HOFF.3HUN  / 

JVChV/2HQu,2HRU,2HU  ,2HV  ,2HW  / 


i19A4 


2007 

200a 

2009 

2010 


20  U 


2012 

2013 

2014 
201b 


2016 


2017 


2016 

2022 


2023 


2024 


2025 


2041 

20A2 

2033 


2034 


FORMAT ( 1014) 

FORMAT (7(A4,A1,3X) ) 

F 0 R M A T ( 14 « 19 A4  ) 

FORMAT(2I4,3F0.3,2I4) 

FORMAT (6F8.3, 14 ) 

FORMAT ( 1oF0,3) 

FORMAT ( A4 i A1 ,515 ) 

FORMAT( 1H0.A4, A1.26H  CANNOT  BE  FOUND  OR  OPENED) 

FORMAT (44H0jUB  ABCRT  - WAKE  PROGRAM  IS  IN  RESTART  MODE) 
iFOKrtAmSHiAKAP  STRATIFIED  SUBMARINE  WAKE  PROGRAM 

FORMAT ( 12H0HUN  RESTART , 11X , 4 A2  ) 

FORMAT ( 10H0RUN  START , 13X , 4A2 ) 

FORMAT (20H0RUN  SPECIAL  REST ART , 3X , 4 A2 ) 

FURMAT(/3X,9HMAX  STEPS, 3X , 10HSTART  STEP , 4X , 5HMIN  X.7X.5HMAX  y 

1 4X,10H1NITIAL  DX,4X,7HMAX  HRS,4X,0HLAM  FUG/  X. 

2 I0»H2,4X,3E12. 4,I8,H2) 

lF?^iT;iru.Iu^E’2X’lOHL01' Y x 

F0RMAT(12X,7E12,4) 

FORMAT (14H0SCALE  F AC  TORS , 4X , 5 < A2 , 1UX ) ) 

FORMAT (15H0WEIGHT  FACTORS , 3X , 5 ( A2 , 1 nX ) ) 

F0RMAT(/3X,9HN0ISt  MIN , 2X , 10HMAX  CHANGE , 2X , i OHEDGE  TOLER. 3X 
1 8HMIN  CURV,4X,8HMAX  CURV/5E12.4)  * lunEUoE  TOLER, 3X, 

1^1  OHFULL^PRiNT^ax  5ohI^L'‘,X’9HINT  PRINT^*!shIM  DISK.!;,’ 

* 3Sa»HAU«1S:Ss;;«:^c?^K;f2x;n;1T2°i4L  print^'8h!*  m,, 

FORMAT(12HOX  OUTPUT  AT,10E12,4) 

FURMAT(/6X,2HKE,10X,lHG,llX,2HPR,8Xt7HSMALL  A,7X,lHA.llX  1HB 
1 10X,4HBETA,9X,1HC,11X,1HS/9E12.4)  »7X,1HA,11X,1HB, 

FORMAT dSHOSOLu-nON  FILES, 4X,7HW0RKING,4X,5HGAMMA.5X.8HGL0Bal 
1 4X , 4hFL0T , 6X , &HCOMMON/20X , 7 ( A 4 , Al , 5X ) ) * 

FORMAT (18H0REFLEC I ION  VECTOR , 2X , 1416 ) 

FORMAT (12H0EDGE  VALUES, 6X,5(A2,10X) , 2HSL , IqX , 2HP  ) 

FORMAT  (14HOCONTROL  F LAGS  , 5X  , 9HT YPE  STAT  , 4X  , 8HSTART-UP  sy 

1 5H3REAK, 6X ,6HoU  1NI,6X,7HBU  FULL,5X,8HBU  TOTa*  4X  8hl T£Ra TES 

2 4X.16H0ISK  SAVE  VEC T OR/12X , 7 1 12 , 9X , 712)  T0T AL ’ 4 X ’ 8hI T£RA F£S • 
FORMAT (17H0STAKT  CONDITIONS , 3X , 4HNAME , 6X , 5H VALUE , 5X , 6H0PT ION , 


0 
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1 4X,6HLAM6DA,3X,8HPKESSURE,3X,6HOUTPUT/20X,A4*Al»l9,4ll0> 
203b  FORMAT (12HOORDZ  CHANGE ,2El2 .4  ) 

204Q  FORMAT(/3X,8HO  CUTOFF, 3X, 10 HD IVG  PCENT , 3X , 9HDI VG  FACT,4X, 

1 6HP  N0RM,5X,6hS  CUTOFF ,4X,9HCJT0FF  MU,3X,8HX  FACTOR, 5X, 

2 6hX  ZEK0,5X,8hY  UUTPUT,4X,8HZ  OUTPUT/10E12.4) 

C 

C RUN  INITIALIZATION 
C 


C 

C 

C 


5 

10 

20 


IOLAYrl 


CALL 
CALL 
CALL 
CALL 
NC0MT=419 
NWwZF=2 


I IBFK ( 2 « BUFR ) 
L£TlI(LTID,BUFK) 
IIBFR<2«BUFS) 
LLTLl(FT10,bUF!>) 


NWR=24 

NMR=lfl 

NM0VE=12 

NMAI=16 

NINU=6 

N0uT=5 

MVAR=5 

NWV£C=o 

NVART=7 

IROWTrl 

IR0WG=6 

IMAPVClIsl 

IMAPV(2)r7 

IMAPV(3)=l3 

IMAPV(4)=19 

CALL  GOATE(IDATC) 


open  common  save  disk  pile  and  read  in  contents 


L=1 

CAUL  LETLU(C0MID,2,BUFR,ET1D,C0MNU)  ,-l,2,N) 

IF  (N)  5,20*5 

CALL  SFVMV(C0MN,FMME,2) 

WRITE  (NOUT ,2000 ) FNAME (L ) ,FNAME (L+l ) 

CALL  EXIT 
NRX=1 

CALL  PBFDR(C0M10,NRX,NC0MT ,NSTST) 


C 

C READ  initial  input  card 

c 

READ  ( NlNU,  1 0 0 0 ) INPLG , N , JOBE 
LSTFLsl 

inflg=inflg-i 
IF  (INFLG)  110,115,120 
110  LSTFL=-1 

NSTSTzNPTS 
GU  TO  1211 

lib  IF  (LENJF)  118,116,118 

lib  IF  (N)  118,117,118 

117  WRITE  (NOUT, 2002) 

CALL  EXIT 

118  LENDFsO 

c 

C READ  REMAINING  INPUT  CARDS 


! 


< 
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i 


c 

120 


J I 


READ 

( NINU , i002 

REAU 

(NINU, 1003 

REAu 

(NINU, 1000 

read 

(NINU, 1000 

read 

(NINU, 1005 

READ 

(NINU, 1005 

Rt.AU 

( NINU, l 005 

Rlau 

(NINU, 1005 

read 

(NINU, 1005 

reau 

(NINU , luOo 

read 

(NINU, 1005 

read 

(NINU, i000 

READ 

(NINU, 1003 

read 

(NINU, 1005 

REAU 

(NINU, 1004 

REAU 

(NJ.NU  ,100d 

1 TOUT 

, ZOUT 

NRUNI , CMNT 

RSlMX,NSTSl ,XMIN,XMAX»DELX,MXhRS,LaMIN 

FXKT,LYLFF,MXRZ,LZLFF 

NRFV 

VSCAV 

VWTFV 

^E«UV 

i.PSN,EPSX,EPSS,ECMN»ECMX 

DXMAX,OXM1N,DXFMX,bUFAC»FCUR,OFFMN,DFFMX,DFR'i!X 

MOLP,NIUPP,NFOLP,NTOPP,NTOLP,NSOLP,NIAAF,LOjT 

XOUTV 

NSTAT,NSlaC,NUFF,LiaUF,LFtjUF,LTBUF,NSTpR,IDSY 
RE*G,PR,aS, a,8,BETa,C,S 
DC*U,Ki5,URDZ,SCALE*PCRIT,SCALM«SP0RS,CP0K 
^>1  • S2 , S5,b6  , S7  , S0 , LSc AL 
«CUT ,DIVP,DIVF,PNORM,PCUT,CMU,XFaCT,XZERO, 


C 
C 
C 

1211 


Rfc.Au  CMlNUflOOl)  FNAME 
OPEN  REMAINING  BASIC  DISK  FILLS 


•1*2*NWR,N) 


L=1 

CALL  LETLU(SlN10,2,BUFS,FUD,FNAME(1) ,< 

IF  ( N ) 10  ♦ 1? 12  » 10 
1212  L=3 

CALL  LETLU(bAMlD,2,bUFK,ETlD,Fi\lAME(3>  ,-l,2*NMAT,N) 
IF  (N)  10*1213*10 


u 


C 

c 
c 

1213 

1214 


open  GLOBAL  disk  FILE  AND  PROFILE  DISK  FILE 
121b, 1215*121*+ 


IF  (NIOPP) 

L = 5 

CALL  LET;.U(bLOlO,2,BUFK,ETID,FNAME(5)  ,-l,2,N) 
IF  (N)  10*1216,10 

1215  CALL  SFVFL ( BLANK ,FNAM£ ( 5 ) f 2 ) 

121b  IF  (NTOPP)  1216, 121«, 1217 
1217  L=7 

CALL  LETLU(PLTlO,2,bUFR,tTlD,FNAME(  7)  ,-i,2,(\j) 
IF  (N)  10,122,10 
CALL  SFVFL(BLANK,FNAME(7) ,2) 

IF  (INFLG)  144,123,122*+ 


1216 
122 
C 

C SPECIAL  RESTART  WITH  CURRENT  WORKING  FILE 
C 

1224  NSTSTrNPTS 
LSTFLr-1 
GO  TO  142 
C 

C PROFILE  INPUT  INI  1 1 ALUAT ION 
C 

NSTSTzNSTSI 
NPTS=NSTST 

IF  (NRUNI)  125,125,124 

nkun=nruni-i 
nkun=nrun+i 

IKOWR=2 
ZSCALrSC ALE 

yscal=scale 


(NE*  COMMON  DATA) 


12A 


124 

125 
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ZEROV(6)=2.0*SCALE 

X=XMIN 

IF  (OELX)  140 , 130 « 129 

129  DX=OELX 

130  READ  ( N I NU *1020 ) F ILEN . LZFAF ,LYFAF ,LZMAX , LPRFL, I YpSN 
IF  (LYFAF)  134,141,141 

131  CALL  L£TLU(F1LI0,2,BUFR,ET1D,FILEN(1) ,-l,2,N) 

IF  ( N ) 132,140,132 

132  WRITE  ( IMOUT  , 2000  ) F1LEN 
CALL  EXIT 

134  CALL  SFVMV(FNAM£(7)  ,FIL£l\l,2) 

C 

c initialize  fixed  parameters 
c 

140  NSS=1 

LlAAFrl 

DEPSTrO.O 

DEPSIsO.O 

OXSAVsO.O 

142  DU  143  1=1, 6, 2 
KLOK ( 1 ) =0 

143  CONTINUE 

ldruz=o 

IF  (ORUZ)  144,1435,144 
1435  LURDZrl 

READ  ( NlNU , 100b ) URUZL,DRDZH 

144  LSOLF=0 

DU  145  1=1,6 
LFLA6 ( 1 ) =0 
lHb  CUN  I 1 NUL 

JERRzJERRX 
FMU=1 • O/Rt 
XK=PR/Rt 

BBETA=1.0+2.0*B*BETA 
HBETA=1.5*BBETa 
BriS  = l , 0+BBETA/d/S 

CVV=( 1.0-2. 0*8* (1. O-BETA) )/3.0/BBETa 
C 

c reflection  property  CONDITION  CONSTRUCTION 
c 

DO  155  L=1 , 2 
DO  151  I=l,NwVtC 
IRFV( I,L)=NRFV( I,L) 

151  CONTINUE 

DO  153  1=2,3 
DO  152  Krl,2 
N=2*(3-i  )+K 

j;RFV(N,L+2)=IAbS(l-NRFV(I,L)-IABS(L-K)  ) 

152  CONTINUE 
155  CONTINUE 

DO  154  N=5,NWVEC 
IRFV ( N, L+2 ) =NRF V ( 7 , L ) 

154  CONTINUE 
15b  CONTINUE 

CALL  OATSW ( 14 , 1SW ) 

SO  TO  (160,200},isw 
160  WRITE  ( NOUT , 1000 ) IRFV 
C 

C OUTPUT  RUN  PARAMETERS) 
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c 

200  WHITE  (NOUT «200b ) CMNT,NRUN 
IF  (INFLG)  201,202,203 

201  WRITE  (NOUT ,2007)  IUATE 
GO  TO  204 

2 02  WRITE  ( NOUT , 2006 ) IUATE 
GO  TO  204 

203  WRITE  (NOUT, 2009)  IUATE 

204  WRITE  (NOUT, 2010)  NSTMX , NSTST , XMIN , XMAX ,DX ,MXHRS ,lAM1N 
WRITE  (NOUT, 2011)  MXRY , LYLFF , MXRZ , L2LFF 

WRITE  (NOUT, 2031)  N«FV 
WRITE  (NOUT, 2013)  JVCHV 
WRITE  (NOUT, 2012)  VSCAV 
WRITE  ( NOUT , 2014 ) JVCHV 
WRITE  ( NOUT , 20 i2 ) VWTFV 
WRITE  (NOUT, 2032)  JVCHV 
WRITE  ( NOUT ,2012 ) ZERO V 

WRITE  (NOUT, 2015)  E^SN , EPSX , ePSS , ECMN, ECMX 

WRITE  (NOUT  ,20lfa)  0XMAX,UXMIN,DXFi4X,BUFAC,  FCUR  ,DFFMN,DFFMX,DFRNX 
WRITE  (NOUT, 2017)  NIOLP ,NI0PP ,NFOLP , NTOPP , NTOLP ,NSOlP , Nl AAF , LOUT 
WRITE  (NOUT, 2016)  XUUTV 

WRITE  (NOUT, 2033)  NSTaT , NST BC , NUF F , Li BUF , LFbUF , LTbUF , NSTPK , I DSV 
WRITE  (NOUT, 2022)  RE , G , PR , AS , A , B , BET A , C , S 

WRITE  (NOUT, 2023)  OC , U, RIS , DRDZ , SCALE , PCR1 T .SCALM, SPURS, CPOR 
WHITE  (NOUT, 2024)  Si , S2 , S5 , S6 , S7 , S6 , HSC AL ( LSCAL+1 ) 

WRITE  (NOUT, 2040)  QCUT ,QI VP,DI VF , PNORM , PCUT , CMU , XFACT ,XZERO, 

1 YOUT , ZOuT 

WRITE  (NOUT, 20^5)  FNAME,COMN 
IF  (INFLG)  210 ,205,210 

20b  WRITE  (NOUT, 2034)  F1EEN,LZFAF,LYFAF,LZMAX,LPRFL, IYPSN 
210  IF  (LDRDZ)  215,220,215 
2lb  WRITE  (NOUT, 2035)  DHUZL , DRDZH 
C 

C STmRT  run  execution 

c 

220  LMLFL=-1 

NPTSN=NPTS 

RETURN 

END 

CART  10  0105  OB  AOOR  5570  DB  CNT  0224 
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WAKIY.S<0105> 

**WAK 1 Y - STRATIFIED  SUBMARINE  WAKE,  INSERT  A Y POINT 
SUBROUTINE  WAKIY(OYM, JTT, JTP,IY) 

c 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  INSERTS 
Cay  point  between  two  existing  ROW  POINTS 
c 

♦CUPY  (CMWAK) 

c 

RATYr(YNEWV(IY)-YNEWV(  IY-1)  ) / < YNEWV < I Y+l ) - YNEWV < I Y-l ) ) 
DO  102  1=1,2 

CALL  WAKMP< I ,NPOS»JYM,XZMV,i) 

CALL  WAKMP(1,NP0S,JYH,XZPV*1) 

DO  101  U=l,NWVcC 

XV  ( J)rXZMV( J)+KATY*(XZPV( j)-XZMV( J) ) 

101  CON  rlNUE 

CALL  L'AKMP  ( I , NPCS  • JY  T , X V , 2 ) 

10k  CONTINUE 
RETURN 
END 

cart  id  o io 5 db  ado:,  i-tco  ob  cnt  0026 
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WAKIZ.S  ( 0105) 

**WAK1Z  - STRATIFIED  SUBMARINE  WAKE  , INSERT  A Z ROW 
SUBROUTINE  WAKIZ < NRM ,NKT , NKP « IZ ) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  INSERTS 
C A Z ROW  BETWEEN  TwO  EXISTING  ROWS 
C 

♦COPY  (CMWAK) 

C 

C FILL  MINUS  aND  PLUS  ROWS  IN  BUFFER 
C 

NP0S=NRM+1 

CALL  WAKRR ( NRM, ZMV ) 

NP0S=NRP-1 

CALL  WAKRR ( NRP , ZPV ) 

Z=ZNEWV(IZ) 

RATZ=(Z-ZM)/(ZP-ZM) 

IF  (IYPSM^IYPSP)  10.10.20 
10  IYPS=MMAX(IYPSM.IYPSP) 

GO  TO  30 

20  IYPS=MMIN( IYPSM.IYPSP) 

30  I YPE=MMAX(IYPEM, IYPEP) 

IYPSX=IYPS+LYLFF 

C 

C INTERPOLATE  FOR  NRT  ROW 
C 

NPOS=NRT 

DO  114  JY=IYPSX,IYPt 

00  112  1=1.2 

CALL  WAKMP(l.NRT-lfJYfXMYV.l) 

CALL  WAKMP(I,NRT  + 1.*JY  » XP Y V , 1 ) 

DO  HI  Jrl.NWVEC 

XV/(d)=XMYV(J)+KATZ*(XPYV(J)-XMYV(J)) 

111  CONTINUE 

CALL  WAKMPd.NRT.dY. XV, 2) 

112  CONTINUE 
114  CONTINUE 

CALL  WAKWR(NRT,ZV) 

RETURN 
END 

cart  id  0105  db  ador  2Lfo 


DB  CNT  004E 
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WAKLL.S(0105) 

♦♦WAKLL  - STKAT IFIED  SUBMARINE  WAKE»  LINER  LOCATION  CALCULATION 
SUBROUTINE  WAKLL(1EK*LFL) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  DETERMINES  WHETHER  THE  LINER 
C HAS  BEEN  KLaCHEO  AND  SLTS  THE  COEFFICIENT  APPROPRIATELY 
C 

♦COPY  (CMWAK) 

C 

TEM=0.0 

RH=Y*Y/YSCAL/YSCAL+Z*Z/ZSCAL/ZSCAL 
IF  (KR-SpORS)  100«1UU«10 
lo  TEMM=cPOR*(RK-SPOKS) 

GO  TO  ( 20  « 30 ) « LFL 

20  TEM=0.5*(£XP(TtMM)+LXP(-TEMM) )-1.0 

RETURN 

30  TEM=Y*(XV(4)-2ER0V(H) ) /YSCAL/YSC AL+Z* < XV < b) -ZEROV < 5) ) /ZSCAL/ZSLAL 

TlM=cPOR*TEM*(lXP(TLMM)-EXP(-TEMM) ) 

100  return 

END 

cart  id  0105  ob  aodk  2E20  ob  cnt  0026 
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WaKMG.S( 0105) 

♦♦WAKMG  - STRATIFIED  SUBMARINE  WAKE,  MOVE  GAMMA  ROW 
SUBROUTINE  WAKMG(NRX,ZPOS,LFL) 

C 

C THIS  SUbROUTlNE  IN  THE  WAKE  PROGRAM  MOVES  THE  GAMMA  MATRIX  ROW 
C BETWEEN  THE  ROW  BUFFtR  AND  ulSK  FILE 
C 

DIMENSION  ZPOS ( t ) 

♦COPY  ( CMWAK ) 

C 

DATA  jERRX/2HMb/ 

C 

IF  (NRX)  100,100,10 
10  IF  (NRX-MXRZ)  20,20,100 

20  CALL  SFVMV(2P0S,ZA,NWWZF) 

IYPSX=IVECA+LYLFF 
LL=NMaT*( IVEC8-IYPSX+1) 

NKXX=(NRX-:  ♦MXRY+IYPSX 
GO  TO  (30,40),LFL 

30  CALL  PBFDK(GAM1D,NRXX,LL,R0WG(1,IYPSX) ) 

RETURN 

40  CALL  PBF0W(GAM1D,NRXX,lL,R0WG(1,IYPSX) ) 

RETURN 

100  JERR=JEKRX 

RETURN 
END 

cart  id  0105  ob  aour  2060  db  cnt  0032 
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WAKMP.S(0105) 

♦♦wakhp  - stratifieo  submarine  wake*  move  a point 

SUBROUTINE  WAKMP( 1MAP,NKZ»NRY,  VEC.LpL) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  MOVES  SELECTED  POlNl 
C INFORMATION  BETWEEN  THE  ROW  OR  GAMMA  BUFFER  AND  VEC 
C 

DIMENSION  VEC ( 2 ) 

♦COPT  (CMWAK) 

C 

DATA  JERRX/2HMP/ 

C 

J=NR2-NP0S+2 
IF  (J)  200.200,10 
10  IF  (J-3)  20,20,200 

20  IF  ( IMAP-5 ) 60 , 20U , 30 

30  GO  TO  (4o,5o)»LFL 

10  CALL  SFVMV(ROWG(l,NKT) ,VEC.NMAT) 

return 

50  CALL  SFVMV(VEC,R0W6(1.NRY) ,NMAT) 

return 

60  NRXXrlMAPV(lMAP) 

IF  (IMAP-3)  90,70.00 
70  IF  (MOOD)  90,60,200 

80  NRXX=NRXX-NWVEC 

90  GO  TO  (100,110) ,LFL 

100  CALL  SFVMV(ROWd(NHXX,NRY, J) ,VEC,NWVEC) 

RETURN 

HO  CALL  SFVMV(VEC,R0W9(NRXX,NRY,  J)  ,NWVEC) 

RETURN 

200  JERRs JERRX 

RETURN 
ENU 

CART  ID  0105  Ob  AOUR  2UA0  OB  LNT  0040 
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WAKMR.S(0105) 

**WaKmR  - STRATIFIED  SUBMARINE  WAKE,  MOVE  A BUFFER  ROW 
SUBROUTINE  WAKMR(MRF *NRT) 

THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  MOVES  THE  CONTENTS  OF  THE 
ROW  BUFFER  FROM  ROW  NRF  TO  ROW  NRT 

♦COPY  (CMWaK) 

C 

00  100  I=1,MXRY 

CALL  SFVMV(ROWb(l,I»  •“RF  ) ,RUWB(1,I,NRT)  » NMR ) 

100  CONTINUE 
RETURN 
END 

CaRT  ID  0 1O5  OB  AODR  2040  DB  CNT  OOlA 
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WaKMV , S ( 0105 ) 

♦#W aKMV  - STRATIFIED  SUBMARINE  WAKE , SUPEREQUILIBRIUM  MATRIX  VALUES 

subroutine  wakmv(Xxv«ttv» n 

c 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  COMPUTES  THE  SUPeREQU^LIBRIUM 

c matrix  values  needed  for  the  implicit  upsweep  calculation 
c 

DIMENSION  XXV(fa) tlTV(6) »T(10) 

♦COPT  (CMRAK) 

C 

IF  ( XX V ( 1 ) ) 100,100*101 

100  CALL  SFVFL(0.0*T«10) 

CVS=0.0 

CwSsC.O 

RETURN 

101  0=SQKT(XXV(1) ) 

FrLM=TTV(4)+DR0X 
CALL  WAKCL(XXV,FTLM) 

IF  (LAMIN)  100*102,100 

102  BOL=BBETA^Q/SCALE 
C1=A^BQL*Q/SCALE-G*F 1EM 
C2=BBS^G*FTEM-Cl 

TEMsTTV ( 1 ) *TTV ( 1 ) *SCALE/Q/BQL-BBET  A^G^FTEM/C2 

1 -TTV(2)«TTV(2)*BUElA*(G*FTEM*(l«0-A/B/S)+A*A*XXV(l) 

2 /SCALE/SCALE)/C1/C2 
CVS=CVV+(B-CVV^TEM)/(HBETA+TEM) 

CALL  WAKTC(CVS,1.0,IVECA) 

C VWrCVS/Cl 
C4=CVW*G*TTV(3) 

T(7)=-dQL*CVW^XXV(l) 

C3=CVS*UR0Z+C4^BbS*FTV(3) 

T (4)=CVS*XXV(1)*BUL/C2 
T(3)=C3*BUL/C2 

CWS=CVS-2»0^G*(C3+CVS*TTV(4))/C2 
CALL  WAKTC(CWS«1.0-CVS,IVLCB) 
DUR=(FTEM*CWS*XXV<l)-*-C4*TTVl3)*XXV(l) 

1 -BQL«(T(3)«XXV(1)*1 ( 4 ) ♦TTV ( 4 ) ) ) /Cl 
T { 2 ) =- ( G*DUR+C WS*XX V ( 1 ) )/bUL 
T(6)=-CVS*XXV(l)/tiQL 
T(5)=-C4^TTV(2)/60L 

CURs(FTEM*C4*TTVil)*TTV(3)*TTV(l)*CVS*(l*0 
1 +BQL*BQL/Ci) )/Cl 
T(l)r-(G*CUR+C4*TTV(1) J/BUL 
T(8)=CVS 
T(9)=CWS 
T ( 10 ) =C4 

return 

end 

CaRT  10  0105  DB  AOOR  41C0  DB  CNT  OObE 
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WAKMY.S(0105) 

**wakhy  - stratified  submarine  wake,  y implicit  MATRIX  COEFFICIENTS 
subroutine  wakmy 
c 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  COMPUTES  THE  X,Y,Z  AND  D MATRIX 

c coefficients  for  the  y implicit  upsweep 
c 

(DIMENSION  T ( 10 , 5 ) 

*CUPY  (CMWAK) 

EQUIVALENCE  ( T ( 1 , 1 ) * ZNEWV ( 1 ) ) 

C 

FPRZ( aRGMi ARG« A«GH)=FZM*( aRGM+FZS*aRG-FZR*aRGP) 

FPKY ( aRGM»  ARGi aRGH ) =F YM* ( ARGM+F YS* ARG-F YR*aRGP ) 

SpRZ( ARGMf ARGi AKGH)=SZM*AKGM+SZ*ARG+SZP*ARGP 
SPRY(  aRGM,ARG,aRGH)=SYM*ARGM-*-SY*ARG  + SYP*ARGP 
c 

c compute  spacing  factors 
c 

DZM=Z-ZM 

DZP =ZP-Z 

OZT=DZM+DZP 

FZM=-dZP/UZM/OZT 

FZp=OZM/OZP/DZT 

FZ=-FZM-FZP 

TEM=L>ZM/DZP 

fzr=tem*tem 

FZS=FZR-1.0 

SZM=2.0/DZM/OZT 

SZP=2.0/UZP/UZT 

SZ=-SZM-SZP 

c 

dym=y-ym 

DYP=YP-Y 

DYT=UYM*DYP 

FYM=-DYP/UYm/DYT 

FYP=UYM/OYP/UYT 

fy=-fym-fyp 

TEM=OYM/DYP 

FYR=TeM*TEM 

FYS=FYK-l.O 

SYM=2.0/DYM/UYT 

SYP=2.0/DYP/UYT 

sy=-sym-syp 


C 


C 

C 

C 


CALL  SFVFL(G.OiXMATiNWVEC) 
CALL  SFVFLQ.OiYMATiNWVEU 
call  SFVFL(0.0,ZMAT,NMOVE) 

COMPUTE  MULTIPLICATIVE  FACTORS 

CALL  WAKOG(ZM) 

CALL  WAKMV(XMYViTP'YViT(lib)  ) 
CALL  WAKDG(ZP) 

CALL  WAKMV(XPYV,TRYV»T(1,4) ) 
CALL  WAKOG(Z) 

CALL  WAKMV(XZMv,TZMV»T<1»3) ) 
CALL  WAKMV(XZPV,TZPV,T(1,2U 
call  WAKMVtXV.TV.T <1i1)  ) 
tmu=fmu 
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0 


u 
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c 

c 

c 


c 

c 

c 


IF  (SROKS)  40.40,20 
TEM=Y*Y/YSCAL/TSCAL+Z*Z/ZSCAL/ZSCAL 
IF  (TE.M-PCUT)  40,40,30 
TMU=CMU 

CL=2.0*(AS*TMU/SCALt+B*G) /SCALE 
COMPUTE  derivative  factors 


OAX=DXI* ( U+XV ( 5) ) 

CALL  WAKLL(TEMK.l) 
AOM=FYM*XZMV(4) 
AO=FY*XV (4) 
ADP=FYP*XZPV(4) 

SL=SC ALE/3. 0 
FGY=FPKY(SORT(XZN,V(l)  ) 
FUZ=FPKZ(SQRT(XMYV(1) ) 
FKZ=FPR2(XMYV(1) ,XV(1) 
FKZ=FpRZ(XMYV(2) ,XV(2) 
FUY=FPRY(XZMV(3) ,XV(4) 
FUZ=FPRZ(XMYV(a) ,XV(3) 
FVYrFPRY(XZMV(4) ,XV(4) 
FVZ=FPRZ(XMYV(4) ,Xv(4) 
FWY=FpRY(XZMV(5) ,XV(5) 
FWZ=FPRZ(XMYVO)  ,XV(5) 
FLY=FPRY(XZMV(fa) ,Xy(6) 
FLZsFPRZ ( XMYV ( fa ) , aV ( b ) 


Q, SORT ( XZPV ( 1 ) ) ) 
Q , SGRT (XPYV ( 1 ) ) ) 
XPYV(l) ) 

XPY V ( 2 ) ) 

XZPV ( 3 ) ) 

XPYV ( 3 ) ) 

XZPV ( 4 ) ) 

XPYV ( 4 ) ) 

XZPV (5) ) 

XPYV (51) 

XZPV (6) ) 


t , 


C 

c 

c 

bO 


XPYV ( 6 ) ) 

SPLY=3.0*OC*YSLAL*CVS*Q+TMU 
FPL  Y = 3 • 0*DC* YSCAL*C VS*FU  Y 
SPLZ  = 3.0*DC*SCAL*tWS*Q+TliO 
FPLZ=3.0*DC*SCAL*LWS*FQZ 

PROM=FUY*T(5,l)+FUZ*r(l,l)+FVY*T(8,l)  + (FWY4-FVZ)*T(10,l  )+FWZ*T(^.l) 
PK0P=FUY*FUY*T(6,1)+FUZ*FUZ*T(2,1)-Q*SCALE»(FWY+FVZ)**2/3.C 

COMPUTE  MATRIX  ELEMENT^  FOR  UU 

TEM0=SPLZ*SPRZ(XMYV(1)  ,XV(1) ,XPYV(1) )+FPLZ*FKZ 
TEMC=FPRZ(XMYV(1)*XMYV(5) ,XV(1>*XV(5) ,XPYV(1)*XPYv(5) ) 

XMAT(1 )=-SPLY*SYM-FPLY*FYM+ADM 

YMAT(D=DXX-SPLY*SY-FPLY*FY  + AD+CL+2.0*(PROM+G*T(3,1)  ) 

ZMAT ( 1 )=-SPLY*SYP-FPLY*FYP+ADP 

DVEC(1)=DXX*XV(I)-2.0*(PROP  + G*FKZ*T(4,1)  h-temd-temc 
COMPUTE  MATRIX  ELEMENTS  FOR  RHO 
TEMM=T(7,1)-XK 

TEMP=FPRY(T(7,3) ,T (7,1) ,1  (7,2) ) 

TEMD=( T(4,1)-XK)*SPKZ(XMYV(2) ,XV(2) , XPYV (2) ) 

1 +FPRZ ( T (4,5) »T (4,1) ,T (4,4) )*FRZ 
TEMQ=FPRZ(T(3,b)*XMYV(l) , T ( 3 , 1 ) *XV ( 1 ) , T ( 3, 4 ) *XPYV ( 1 ) ) 
TEMC=FPRZ(XMYV(2)*XMYV(5) ,XV(2)*XV(5) ,XPYV(2)*XPYv(5) ) 

XMAT ( 2 ) =TEMM*SYM+ I EMP*FYM4ADM 
YMAT (2)=DXX+TEMM*SY+TEMP*FY+AD 
ZMAT(2)=TEMM*SYP+1EMP*FYP+ADP 
DVEC(2)=DXX*XV(2)-TEMD-TLMQ-ORUZ*XV(5)-TEML 
IF  (NSTAT-2)  50,60,50 

COMPUTE  MATRIX  ELEMENT^  FOR  U 

TLMM=T(6,1)-TMU 
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TEMP=FPRY(T<6*3) *T(6«1),T(6,2) ) 

TEMD=(T(2«1)-TMU ) *SPKZ ( XMY  V ( 3 ) « XV ( 3 ) « XPYV ( 3 ) ) 

1 +FPRZ(T<2,5) *T(2«1) «T<2*4) )*FUZ 

TEMQ=FPKY(T(5,3)*XZnV(l) ,T(5*1)*XV<1) «T(5«2)*XZPV<1) ) 

1 +FPRZ(T(1.5)*XMYV(l).T(l,l)*XV(l)«T(lt4)*XPYV(l) ) 
TtMC*FPRZ(XMYV(3)*XMYV(5)«XV(3)*XV(5).XPYV(3)*XPYV(5) ) 
XMAT(3>=TEMM*SYM+TEnP*FYM+ADM 
YMAT  < 3 ) =OXX+TEMM*SY+T£MP*FY+AD 
2NAT(3)=TEMM*SrP+TEMP*FYP+ADP 
DVEC(3>=DXX*XV<3)-TEM0-7EM8-TEMC 
IF  (NSTAT-2)  70,60,60 
C 

C COMPUTE  MATRIX  ELEMENTS  FOR  V 
C 

60  TEMD=Q*Sl*FPRZ <FPKY<XMMV< 5)  ,XMYV< 5) ,XMPV( 5)  ) ,FWY* 

1 FPRY (XPMV(b) «XPYV(b),XPPV(5) ) ) 

TEmO=fPRY( T(8,3)*XZnV(l)  *T<6, 1 )*XV(1) , T ( 8 , 2 ) *XZPV  ( 1 ) ) 

1 +FPRZ(T  <10«5)*xMYv<l) «T<10«1 )*XV<1)  , T ( 10 , 4 ) *XPYV  ( 1 ) ) 
TLMC=FPRZ(XMYV(4)*XMYV(5) ,XV(4)*XV(5) ,XPYV(4)*XPYV(5) ) 
XMAT(4)=“TMU*SYM+ADM 
YMAT(4)=0XX-TMU*SY+AU+TEMK 
ZMAT(4)=-TMu*SYP+AOP 

DVEC (4)=0XX*XV(4)+TMU*SPRZ(XMYV(4) ,XV(4) ,XPYV(4) ) 

1 +TEM0-FPRY<TZMV<6) ,TV<6) ♦ TZPV ( 6 ) ) -TEMQ+TEMK*ZEROV  ( 4 ) 

2 +Q*Sl*SPRZ(XMYV(4) , XV < 4 ) « XP YV ( 4 ) ) +SL*F8Z* (FWY+FVZ ) -TEMq 
C 

C COMPUTE  MATRIX  ELEMENTS  FOR  W 
C 

TEMD=Q*SL*FPRY <FPKZ<XMMV< 4) , XZMV ( 4 ) , XPMV < 4 ) ) ,FVZ* 

1 FPRZ (XMPV(4) ,XZPV(4) , XPP V ( 4 ) ) ) 

TEMQ=FPRY  (T(10,3)*XZriV(l),T(lO,l)*XV(l)',T(lO,2)*XzPV(l)) 

1 +FPRZ(T<9,b)*XMYV<l) , T < 9 , 1 ) *X V < 1 ) » T < 9* 4 ) *XPYV ( 1 ) ) 
TEMC=FPRZ(XMYV<5)*XMYV<5) ,XV(5)*XV<5).XPYV(5)*XPYv<5) ) 

XMAT<  5)=-TMU*SYM+ADM-SL*<8»SYM+FQY*FYM) 

YMaT<5)=DXX-TMU*SY  + A0+TEMK-SL*  <Q*SY+F8Y*FY) 
ZMAT(5)=-TMU*SYP+A0P-SL*<8*SYP+FQY»FYP) 
DVEC(5)=0XX*XV(5)+TMU*SPRZ(XMYV<5' ,XV<5) ,XPYV<5) ) 

1 +TEMD-FPRZ(TMYV<6) »TV<6>  ,TPYV<6) > -TEMQ+TEMK*ZEROV < 5 > 

2 +SL*F8Y*FVZ-G*XV<2)-TEMC 
C 

C COMPUTE  MATRIX  ELEMENTS  FOR  SLALE 
C 

70  IF  <XV(1)-8cUT*FMAXV(i) ) 60,80,90 
60  DVEC(6)=2.0*SCALE 
GO  TO  100 

90  TEMD=SPLZ*SPRZ(XMTV<6) ,XV(6) ,XPYV<6) )+FPLZ*FLZ 

TEMC  = FPRZ(XMYV(6)*XMYV(b)  ,XV(6)'*XV<5)  »XPYV(6)*XPYV<5)  ) 

XMAT ( 6 ) =XMAT ( 1 ) 

YMAT(6)=DXX-SPLY*SY-FPLY*FY+AD 

1 +Sl*(HROM+PROP/XV(l)  )-*>S2*LL/2.0  + S5*G*(T(3,l)+m,l)*FRZ/XV<l)  ) 

2 ♦S7*xV(6)*<F8r*F8Y+»’UZ*FUZ)/8  + S8*<FQY*FLY+F8Z*FLZ) 
ZMAT(6)=ZMAT(1) 

DVEC<6>=DXX*XV(6)+TtM0-S6*8*<FLY*FLY+FLZ*FLZ)-TEMc 

C 

C LOWER  BOUNDARY  CONDITION  CHECK 

C 

100  IF  (JY-2)  200,110*105 

105  IF  (JY-IYPSN)  300*300,500 

110  IF  (LYLFF)  500,300,300 
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200  DO  220  1 = 1 « NWVEC 
TEM=1.0 

IF  ( IRFV  ( I « 1 ) ) 20t> » 205 1 210 
205  TEMx-l.Q 

210  ZMAT(I)=ZMAT(I)+TtM*XMAT(l) 

220  CONTINUE 
GO  TO  400 

300  00  320  1 = 1 « NWVEC 

OVLC ( I ) = DVEC ( I ) -X*Al ( I )*ZER0V(I > 
320  CONTINUE 

400  CALL  SFVFL(0.C»XHAT»NWVLC) 

500  RETURN 
ENO 

CART  ID  0105  OB  ADOR  45C0  0B  CNT  0l7A 


0 

0 

1] 
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*AKMZ.S<0105) 

♦♦WAKMZ  - STKATIFIEO  SUBMARINE  WAKE,  Z IMPLICIT  MATRIX  COEFFICIENTS 
SUBROUTINE  WAKMZ 
C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  COMPUTES  THE  XtY.Z  AND  0 MATRIX 
C coefficients  FOR  THE  Z implicit  UPSWEEP 
C 

DIMENSION  T ( 10  t 5 ) 

♦COPY  (CMWAK) 

EQUIVALENCE  ( T ( 1 « 1 ) « ZNEWV ( 1 ) ) 

C 

FPRZ( aRGM*ARG.ARGH)=FZM^( ARGM+FZS^ARG-FZR^ARGP) 
FPRY(aKGM«ARG.aRGH)=FYMMaRGM+FYS^ARG-FYR^ARGP) 

SPKZ (ARGM«ARG*ARGH ) =szm^argm+sz* arg+szp^argp 
SPRY(A«GM«ARG.ARGP)=SYM*ARGM+SY^ARG+SYP^ARGP 
c 

C COMPUTE  SPACING  FACTORS 

c 

DZM=Z-ZM 

GZP=ZP-Z 

DZT=DZM+OZP 

FZM=-D^P/UZM/OZT 

FZPzOZM/DZP/QZT 

FZ=-FZM-FZP 

TEM=OZM/OZP 

fzr=tem^tem 

FZS=FZR-1.0 

SZM=2.0/0ZM/DZT 

SZP=2.0/0ZP/0Zr 

S*=-SZM-SZP 

C 

dym=y-ym 

dyp=yp-y 

DYT=OYM+OYP 

fym=-dyp/uym/dyt 

FYP=DYM/DYP/DYT 

fy=-fym-fyp 

TEM=DYM/DYP 

FYRsTEM^TEM 

FYSsFYR-1.0 

SYM=2.0/DYM/DYT 

SYP=2.0/0YP/0YT 

sy=-sym-syp 

C 

CALL  SFVFHO.OiXMATiNWVEC) 

CALL  SFVFL(1.0«YMAT«NWVEC) 

CALL  SFVFL(0.0,ZMAT«NMOVu) 

C 

C COMPUTE  MULTIPLICATIVE  FACTORS 
C 

CALL  WAKOG(ZM) 

CALL  WAKMV(XMYV,TMYV«T(1.5)  ) 

CALL  WAKOG(ZP) 

CALL  WAKMV(XPYV»TPYV«T(1»4) ) 

CALL  WAKOG(Z) 

CALL  WAKMV(XZMV»TZMV«T(l«i) ) 

CALL  WAKMV(XZPV«TZPV«T(i»2) ) 

CALL  WAKMV(XVtTVt  MltlJ) 

TMU=FKU 


20 


1 

D 

0 

U 

u 

u 

LI 

U 

0 

U 

D 

I] 

0 

0 
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30 

40 


C 

C 

C 


C 

C 

C 


C 

C 

C 


C 

C 

C 


IF  (SPORS)  40»40«20 
TEM=Y*Y/YSCAL/YSCAL+Z*Z/ZSCAL/ZSCAL 
IF  (TEK-PCUT)  40  ,40*30 
TMUsCMU 

CL=2.0*(AS*TMU/SCALt+B*Q) /SCALE 


COMPUTE  DERIVATIVE  F AC  TORS 


DXX=DXl*(U+XV(d) ) 

CALL  WAKLL(TEMR«1) 
AUM=FZM*XMYV(5) 
AD=FZ*XV ( 5 ) 
ADP=FZP*XPYV (5) 

SL=SC ALL/3 • 0 
FQY=FPRY(SQRT(XZMV(1) ) 
FUZ=FPR2(SQKT(XMYV(1) ) 
FKY=FPRY(XZMV(1) «XV(l) 
FKY=FpKY(XZMV(2) «XV(2) 
FRZ=FPRZ(XMYV(2) «XV(2) 
FUY=FPRY(XZ|SV(3)  ,XV(3) 
FUZ=FpRZ ( XMY V ( 3 ) ,Xv< 3) 
FVY=FPRY(XZMV(4) «XV(4) 
FvZsFpRZ ( XM Y V ( 4 ) ,Xv(4) 
FWY=FPRY(XZMV(S) »XV(5) 
FWZ=FPRZ(XMYV(3) fXv(5) 
FLY=FPRY(XZMV(6) ,XV(6) 
FlZ=F  pRZ ( XMY V ( b ) tXV(6) 


Q « SORT ( XZPV ( i ) ) ) 
QtSURT(XPYV(l)  ) ) 
XZPV(l)  ) 

XZPV ( 2 ) ) 

XPYV (2)  ) 

XZPV ( 3) ) 

XPYV { 3 ) ) 

XZPV (4 ) ) 

XPYV  < 4 ) ) 

XZPV (5) ) 

XPYV  C 5 ) ) 

XZPV (6 ) ) 

XPYV ( 6 ) ) 


SPLY=3.0*DC*YSCAL*CVS*g+TMU 
FPLY=3.0*DC*YSLAL*CVS*FQY 
SPLZ=3.0*DC*SCAL*LWi>*Q+TMu 
FPLZ=3»  0*OC*SCAl*LWS»*FQZ 

PROM=FUY*T(5»l)+FUZ*Ta,l)+FVY*T(Btl)  + (FWY+FVZ)*T(10tl)+FWZ*T(ytl) 
PROP=FUY*FUY*T(6,l)+FUZ*FUZ*T(2tl)-Q*SCALE*{FWY+FvZ)**2/3.0 


COMPUTE  MATRIX  ELEMENTS  FOR  QU 


TEMD=SPLY*SPRY ( XZMV (1) t XV ( 1 ) iXZPV(l) )+FPLY*fKY 
TEMC=FPRY(XZMV(1)*XZMV(4) ,XV(1)*XV(4) »XZPV(1)*XZPV(4) ) 
XMAT(1)=-SPlZ*SZM-FPL2  FZM*ADM 

YMAT(l)=DXX-SPLZ*i»Z-F»  LZ*FZ+AD+uL+2.0*(PROM+G*T(3, 1)  ) 
ZMAT(1)=-SPLZ*SZP-FPLZ*FZP+AUP 

DVEC(i)=DXX*XV(i)-2.0*(PROP4G*FRZ*T(4»1) >+TEMD-TEMC 
COMPUTE  MATRIX  ELEMENTS  FOR  RHO 
TEMM=T ( 4 , 1 ) -XK 

TEMP=fPRZ(T(4«5) « I (4, 1 ) , T (4,4 ) ) 
TLM0={T(7«1)-XK)*SPKY(XZMV(2)  ,XV(2> ,XZPV(2) ) 

1 +FPRY(T(7»3)«T(7»1) «T(7»2) )*FRY 
TEivI(J=FPRZ(T(3«o)*XMYV(1)  ,7 ( 3 « 1 > *XV ( 1 ) ,T (3 ♦ 4 5 *XPY V ( 1 ) ) 
TEMC=FPRY(XiMV(2)*XZMV(4) ,XV(2)*XV(4> , XZPV ( 2 ) *XZPV ( 4 ) ) 
XMAT(2)  = TEMM*SZM+TEMP*FZIifADM 
YMAT ( 2 ) =OXX+TEMM*SZ+TEMP*FZ+AD 
ZMAT(2)=TEMM*SZP  + 1EMP*FZP4-ADP 
0VEC(2)=0XX*XV(2)-TLMQ-TEM6I-DRDZ*XV(5)-TEMc 
IF  (NSTAT-2)  50,60,50 


compute  matrix  elements  for  u 


1 


J 


50  T£MM=T ( 2 , 1 ) -TMU 

T£MP=FPKZ(T(2.5)  ,T  (2,1)  ,T<2«4) ) 

TEM0=<T<6,1)-TMU)*SPKY(XZMV<3) *XV(3) ,XZPV(3) ) 

1 +FPRY(T(6»3) « T ( 6 * 1 ) * T ( 6 «2 ) ) *FUY 
TEMa=FPKY(T(5.3)*XZMV(l) «T(5»1)*XV(1) ♦T(5.2)*XZPV(1) ) 

1 +FPRZ(T(1«5)*XMYV(1) «T(1«1)*XV{1) *T(1*4)*XPYV<1»  ) 
T£MC=FP«Y(XZMV(3)*X2MV(4)  «XV(3)*XV(4)  ♦XZPV(3)*XZPV(‘f)  ) 
XMAT(3)=TEMM*SZM+TEMP*FZM+ADM 
YMAT(3>=DXX+TEMM*SZ+TEMP*FZ+AD 
ZMAT(3>=TEMM*S*P+IEMH*FZP+ADP 
CVEC(3)=0XX*XV(3)-TEM0-TEMa-fEMC 
IF  (NSTAT-2)  70,60,fa0 
C 

C COMPUTE  MATRIX  ELEMENTS  FOR  V 
C 

60  TLMD=Q*SL*FPRZ(FPKY(XMMV(5) tXMYV(5) #XMPV(5) ) »FWY« 

1 FPRY(XPMV(5) »XPYV(5) ,XPPV<5) ) ) 

TEmQ=FP«Y(T(6»3)*XZMV(1) . T < 8 . 1 ) *XV < 1 ) » T ( 6 , 2 ) *XZPV < 1 ) ) 

1 +FPRZ(T(10.5)*XMYv(1).T(10,1)*XV(1).T(10«4)*XPYV(1)) 
TEMC=FPRY(XZMV(4)*XZMV<4) ,XV(4)*XV<4) ,XZPV(4)*XZPV<4) ) 
XMAT<4)=-TMU*SZM+A0M-SL*<e*SZM+FaZ*FZM) 
YMAT(4)=OXX-TMU*SZ+AD+TEMK-SL*(0*SZ+FaZ*FZ) 
ZMAT(4)=-TMU*SZp+A0P-SL*(a*SZP+FaZ*FZp) 
DVEC(4)=DXX*XV<4)+TMU*SPRY(XZMV<4> ,XV(4).XZPV<4)  ) 

1 +TEMD-FPRY(TZMV(6) ,TV<6) ,TZPV<fa) ) -T£MQ+TEMK*ZeROV < 4) 

2 +sl*foz*fwy-temc 

C 

C COMPUTE  MATRIX  ELEMENTS  FOR  W 
C 

TEMQ=Q*SL*FPRY  <FPKZ<XMMV<  4)  ,XZMV<4>  »XPMV(4)  ) »FVZ» 

I FPRZ(XMPV<4) «XZPV<4) ,XPPV<4)>) 

TEMa=FPRY(T<10,3)*XZMV<l) , T < 10 , 1 > *X V < 1 > . T < 10 . 2 ) *XZPV ( 1 ) ) 

1 ♦FPRZ(T(9,5)*XMYV<1)  «T(9«1)*XV(1 ) » T ( 9 , 4 ) *XPYV ( 1 ) ) 
TEMC=FPRY(XZMV(5)*XZMV<4) «XV<5)*XV(4) »XZPV<5)*XZPV<4) ) 
XMAT<5)=-TMU*SZM+A0M 
YMAT(5)=OXX-TMU*SZ+AD+TEMK 
ZMAT(5)=-TMU*SZP+A0P 

OVEC(5)=OXX*XV(5)+TMu*SPRY(XZMV(5) « X V ( 5 ) « XZPV ( 5 ) ) 

1 +TEMD-FPRZ(TMYK6)  « TV<fa),TPYV(6)  )-TEMQ+TEMK*ZEROV<5> 

2 +a*SL*SPRY(XZMV(b) • XV  < 5 > ,XZPV(5) ) +SL*FQY* < FWY+FVZ > -G*XV< 2 ) -TEMC 
C 

C COMPUTE  MATRIX  ELEMENTS  FOR  SCALE 
C 

70  IF  (XV(1)-QCUT*FMAXV<  l)  ) 80,80*90 

60  DVEC(6)=2.0*SCALE 

GO  TO  100 

90  TEM0=SPl.Y*SPRY(XZMV<6)  , XV < 6 ) , XZPV < fa ) )+FPLT*FLY 

TEMC=FPRY(XZMV(6)*XZMV<4) , XV < 6 ) *XV < 4 > » XZPV < 6 ) *XZPV < 4 ) ) 
XMAT(6)=XMAT<1) 

YMAT(6)=0XX-SPLZ*SZ-FPLZ*FZ+AD 

1 +S1*  ( PROM+PROP/XV (1) )+S2*Cl/2»0+S5*G*(T(3,1)+T(4,1)*FR^/XV(1) ) 

2 ♦S7*xV(6)*<FaY*FyY+PUZ*FU*)/a«-Se*(F8Y*FLY+FQZ*FLZ) 

ZMAT ( 6 ) =ZMAT ( 1 ) 

0VEC(6)=0XX*XV(fe)+TEM0-S6*a*<FLY*FLY+FLZ*FLZ)-TEM' 

C 

C lower  boundary  condition  CHECK 
C 

100  IF  (JZ-1)  200  » 110  « 102 

102  IF  <JY=*IYPSM)  104*500,106 


* I 
1 


A-45 


CaKT 


IF  (JY-1)  500*b00<300 
IF  < JY-IYPEM)  b00»50U»300 
IF  (L2LFF)  500»30U.b00 
DO  220  1 = 1 « NWVEC 
T£M=1  • 0 

IF  ( IRFV ( I » 2 ) ) 205 « 20b» 210 
TEM=-1 # 0 

ZMAT  ( I ) =ZMAT  ( I ) -f  TEM*XMAT  ( I ) 
CONTINUE 
GO  TO  400 
DO  320  I=l«NwVEC 

DVEC(I)=DVEC(I)-XMaUI)*ZEKOV(I  ) 

continue 

CALL  SFVFL(0.0,XMAT»NWVEC) 

RETURN 

ENL 

ID  0105  DB  AODR  4740  DB  CNT  0180 
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WAKOT.S(C105) 

**WAKOT  - STRATIFIED  SUBMARINE  WAKE,  OUTPUT  THE  RESULTS 
SUBROUTINE  WAKOT 
C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  OUTPUTS  THE  RESULTS 
C 

DIMENSION  JVCHV(7) ,0P0S(2) «PBUF(660) ,KLOCK(g) ,HDSV(2,2) 

DIMENSION  TOTAL(7«40,2),OVEC(33)  »DBUF(440)  ,NN(40),«M(40)  ,HNULU) 
♦COPY  (CMWAK) 

EQUIVALENCE  ( TOTAL (1,1,1 ) ,ROWG( 1,1)  ) 

EQUIVALENCE  (PBUF(l)  «ROWB  (1,1,3  ) ) , ( PBUF  ( 44l*)  ,OBUF  < 1 ) ) 

EQUIVALENCE  (NN(1) ,NKNZV(l) ) , (MM(l) ,NRNYV(i) ) 

C 

DATA  JP0S/1HY.1HZ/,UERRX/2H  /,HNUL/4H  ,4HNULL/ 

DATA  JVCHV/2HQQ,2HR0,2HU  ,2HV  ,2HW  ,2HSL,2hP  / 

DATA  HQSV/4H  ,4h  , 4H ( STO , 4HRED > / 

C 

1000  F0RMAT</6X,24HFULL  PROFILE  OUTPUT  FOR  , A2 , IX * 2A4 » IX , A4 ) 

1001  FORMAT ( 1H  ) 

1003  F0RMAT(/7X,bHZ  Y,lOEl2.4) 

1004  F0RMaT(/6X,23HMeSH  PROFILE  OUTPUT  AT  tAl»2*  = ,El2.4// 

1  7X,A1,7(10X,A2)  ) 

1010  FORMAT ( / 6X « 52H TURBULENCE  PROFILE  OUTPUT  AT  Z = ,El2,4// 

1 7X(lHY,10Xt2HUV.10X,2HUW,10X,2HVW,10X,2HUR,lOX,2HVR,10X,2HWR, 

2 10X,2HRR,10X«2HUU,10X,2HVV,10X«2HWW/) 

1015  FORMAT ( /6X , 8 HX  POINTS, 9X,1HX,13X, 2 HDK,10X, 

1 8HY  SPREAD, 7X,8HY  PUlNTS, 7X,8HZ  SPREAD, 7X, 

2 8HZ  POINTS, 9X.3HX/U,11X,3HB  V/I 1 0 , 4X , 3E15. 5 , 

3 I9,6X»E15.5, I9.6X.2E14.5) 

1021  FORMAT ( //bX , 6HST  A 1US«7X«8H2  LAMBDA , 7X , 8HR  LAMBDA, 

1 7X.8HY  LAMBDA, 7X,8HM0MENTUM, 

2 7X.8HP  ENERGY, 7X.8HK  ENERGY « 14X , 12HELAPSED  TIME/ 

3 6X,A4,4X,8E15.5,8X,I2,5H  HRS  ,I2,5H  MIN  ,I2,4H  SeC ) 

1036  F0RMAT(/14X,7( 10X,A2) ) 

1037  FORMAT ( 18H  MAXIMUM  VALUE  ,7E12.4) 

1038  FORMAT  ( 18H  Y LOCATlUN  ,7E12.4) 

1039  FORMAT ( 18H  l LUCATION  ,7E12.4) 

1040  FORMAT ( 18H  MAXIMUM  CHANGE  ,7E12.4) 

1041  FORMAT ( 18H  GLOBAL  MAXIMUM  ,7E12.4) 

1045  FORMAT (/18X,2HUV'10*'2HUW, 10 X,2HVW,1 OX, 2HUR,10X,2hVR,10X,2HWR, 

1 10X,2HRR,10X,2HUU,10X,2HVV,10X,2HWW/12H  MAXIMUMS  ,10E12.4) 
1049  F0RMAT(//6X,6HSTATUS,8X,8HX  POINTS , 9X , 1HX/6X , A4 , I 14 , E19 , 5 ) 

2000  FORMAT ( //19H  BACKUP  NEEDED  FOR  ,A2) 

2001  FORMAT ( //26H  DIVIDE  CHECK  HaS  OCCURRED) 

2002  FORMAT ( //22h  OVERFLOW  HAS  OCCURRED) 

2003  FORMAT (//33H  SOLUTION  HAS  SPREAD  TOO  FAST  IN  ,A1,5X,4I6) 

2005  FORMAT ( //22H  RUN  SUSPENDED  AT  X = , E 12 , 5 , 5X , 11HX  POINTS  = ,13, 

1 5X,6HNGS  = • 1 3 » 5X  , 6HNPS  = ,13) 

2006  FORMAT ( //48h  MaXIMuM  POINTS  EXCEEDED  ON  AUTOPOlNT  ADJUST  IN  ,Ai) 

2007  FORMAT ( //22H  GLOBAL  FILE  IS  FILLED) 

2008  FORMAT ( //20H  PLOT  FILE  IS  FILLED) 

2009  FORMAT ( //4H  WAK,A2,19H  ERROR  HAS  OCCURRED) 

2010  FORMAT ( /20H  TOTAL  ELAPSEU  TIME  ,I2,5H  HRS  ,I2»5H  MIN  ,12i4H  SED 

2011  FORMAT ( //49H  TOTAL  PROFILES  STORED  AT  THIS  X ON  THE  PLOT  FILE) 

2012  FORMAT ( //26H  MAXIMUM  RUN  TIME  EXCEEDED) 

2031  FORMAT  ( HE12  • 4 ) 

2032  FORMAT (El2,4,12X*9E12,4) 

2033  FORMAT (El2*4«24X,8E12*4) 

2034  FURMAT(E12.4i36X,7E12.4) 
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2035  F0RMAT(El2.4»4bX,6E12.4) 

2036  FORMAT(El2.4,60X,5E12.4) 

2037  F0RMAT(El2.4,7zX,4E12.4) 

2036  F0RMAT(El2.4,84X,3E12.4) 

2039  F0RMAT(£l2,4«SSX»2E12,4) 

2040  FORMAT(E12.4,108X,E12.4) 

C 

C CHECK  FOR  ERRORS 
C 

IOLAYxlO 

MOOD=l 

NKSTX=NKST-1 

NREND=nKST+JEND 

NZP=JEN0+2 

NYPSX=NYPS-1 

NYPLX=iMYPE  + 1 

NYP=NYPE-NYPS+3 

LTRNF=l-TRNF  + LOVCF+LUVFP+LPKRF-*-LFCUK 
IF  (LBURF)  105,105,100 
100  WRITE  ( >JOUT  « 2000  ) JVCHV ( L60RF ) 

10b  IF  (JERK-JERRX)  106»107«106 
108  WRITE  ( NOUT « 2009 ) JtRR 

LTRNFrl 

107  IF  (LTRNF)  110,200,110 

110  IF  (LDVCF)  120.130,120 

120  WRITE  (NOUT, 2001) 

130  IF  (LOVFF)  140,150,140 

140  WRITE  (NOUT, 2002) 

150  IF  (LPKRF)  160,190,160 

160  WRITE  (NOUT, 2003)  jPOS-(  LPKRF  ) ,NYPS,NYPE,NRST  , JEND 

190  IF  (LFCUR)  195,200,195 

195  WRITE  (NOUT, 2006)  JPOS(LFtUR) 

200  IF  (LPR-CL)  205,900,900 

C 

C CHECK  FOR  INTERMEDIATE  PRINTOUT  TO  LINE  PRINTER 
C 

205  IF  (LIOLF)  210,300,210 

210  CALL  CLOCK(l,KLOCK) 

Y = YOLDV(NYPE'()-YOLDV(NYPSX) 

Z=ZOLDV(NRENO)-ZOLDV(NRSTX) 

YP=XFACT*XP+XZLRO 

ZP=YP*SURT(G)/3.14lb 

WRITE  (NOUT, 1021)  JST AT , ( EPSSV ( I ) , 1=1 » 3 ) »XM0M» XPE , XKt « 
1 ( KLOCK( I ) , 1=1 ,6,2) 

WRITE  (NOUT, 1015)  NPTSN , XP ,DX , Y ,NYP , Z ,NZP , YP , ZP 
WRITE  (NOUT, 1045)  TURBX 
WRITE  (NOUT, 1036)  JVCHV 
WRITE  (NOUT, 1001) 

WRITE  (NOUT, 1037)  FMAXV 
WRITE  (NOUT, 1038)  YMAXV 
WRITE  (NOUT, 1039)  ZMAXV 
WRITE  (NOUT, 1040)  TMAXV 
WRITE  (NOUT, 1041)  GMAXV 


check  for  intermediate  printout  to  GLOBAL  disk  file 

) IF  (LIOPF)  310,400,310 
) NRX=33 

IF  (NCFKL(GL0I0,NG,IMKX+1)  ) 320,330,330 


r 


I 
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320  LTRNFsl 

WRITE  ( NOUT  « 2007  > 

GO  TO  400 
330  OVEC ( 1 ) =XP 

CALL  SFVHV( YSCAL«OVLC(2) »2) 

CALL  SFVMV(XM0M,0VEL(4) t20) 

CALL  SFVMV(GMtOVEC<24) ,8) 

OV  r.C  ( 32 ) =GM  ( 11 ) 

OVEC ( 33 ) =DEPST 

CALL  PBFDW(GLO10,NGtNRX,0VLC) 

NRX=NG 

CALL  PBFDW(GL01D«NRX«l,-i.O) 

c 

C CHECK  FOR  mesh  PR1NTOU1  TO  LINE  PRINTER 
C 

400  IF  (LFOLF)  404 1 450 ,404 

404  WRITE  (NOUT, 1049)  JSTaT,NPTSN,XP 
DO  409  I=1.NVAKT 

IF  ( F M A X V { I ) ) 409 1 409 ,405 

405  DO  406  JY=NYPSX,NTPLX 

TOTAUIt  JYt  J )=TOTAL(  It  JYtl)/FMAXV(  I) 

40b  CONTINUE 

DO  407  JZ=NRSTX,NKEND 
T0TAL(I.JZ.2)=T0TAL(I,dZ,2)/FMAXV(I ) 

407  CONTINUE 

409  CONTINUE 
ZM=LPSSV(l)/2.0/C 
YM=EPSSV(3)/2.0/C 

WRITE  (NOUT, 1004)  JPOS < 2 ) t^OUT , JPOS ( 1 ) , JVCHV 
WRITE  (NOUT, 1001) 

00  410  JY=NYPSX,NYPtX 
Y=YOLDV( JY)/YM 

WRITE  (NOUT, 2031)  Y t ( TOTAL ( I , JT t 1 > t 1=1 t NVART ) 

410  CONTINUE 

WRITE  (NOUT, 1004)  JPOS ( i ), YOUT , JPOS ( 2 )t  JVCHV 
WRITE  (NOUT, 1001) 

DO  420  JZ=NRSTX,NRENU 
Z=ZOLDV( JZJ/ZM 

WRITE  (NOUT, 2031)  Z t ( TOTAL < I , JZ t 2 ) t 1=1 t NVART ) 

420  CONTINUE 
C 

C CHECK  FOR  TOTaL  TURBULENCE  PRINTOUT  TO  LINE  PRINTER 
C 

450  IF  (LSOIF)  460,500,460 

460  IF  (LAMIN)  500,461,500 

461  WRITE  (NOUT, 1049)  jSTAT ,NP1 SN, XP 
GMAX=1.0E-06*FMAXV(1) 

UO  480  UZ=NRSTX,NR£NO 
IYPS=I YPSV ( JZ ) 

IYPE=IYPEV( JZ) 

IF  (IYP5)  480 1 480  t 462 

462  IYPSX=IYPS4LYLFF 
CALL  WAKMG< JZtZVtl) 

1 YS=0 

DO  468  JY=I YPSX , I YPE 
DO  465  I=ltlO 

IF  (ABS(ROWG( It JY) )-GMAX)  465t46bt466 
465  CONTINUE 
GO  TO  468 


( 
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612 

6121 

6122 
6125 
6i3 
61*2 
61*5 


i IY£=JY 

IF  (IYS)  468,4b7,468 
1 IYSsJY 
s COiMTlNUE 

IF  (IYS)  469*460*469 
? WHITE  ( NOUT  ,1010)  ZOLDV(JZ) 

DO  470  JY=IYS,IY£ 

WHITE  (NOUT, 2031)  YULOV ( JY ) , ( ROWG ( 1 , JY ) , 1=1 , 10 ) 
3 CUN  I If  UE 
3 CONTINUE 

CHECK  FOR  TOTAL  PRINTOUT  TO  PROFILE  DISK  FILE 

3 IF  (LTOPF)  510,600,510 
J mCA  = lMYP 

I VtCBrNZP 
J = 0 

00  515  1=1 ,NVART 
J=J+IDSV( 1 ) 

3 CONTINUE 

NRX=J*IVECA*lVECB+IVECA  + lV£Cb+-> 

IF  (N0FHl(PLT1u.NP,NHX) ) 520*530, 530 
3 LTRNFrl 
LTOPFrO 

PHI TE  (NOUT, 2008) 

GO  TO  oOO 
) ZArXP 

CALL  P0FQW(PLT1u,NP»2,Za) 

CALL  P0FOW(PLTIO,NP,1VECA,YOLDV(NYPSX) ) 

CALL  P0FDW(PLTIO,NP,lVEC0.ZOLDV(NRi>TX)  ) 

WHTTE  (NOUT, 2011) 

GC  TO  601 

CHECK  FOR  TOTAL  PH1NIOUT  TO  LINE  PRINTER 

) IF  (LTOLF)  faO 1 , 70  0 , bU 1 
L IK0WR=2 

IF  (LEURF)  604,b04,6U3 
i IK0WR=4 

> IF  (LTOLF)  608*610,608 

J WHITE  (NOUT, 1049)  JbTAT , NPTSN, XP 
) DO  680  IVAH=1*NVAHT 

CALL  SFVFL(2EKUV(IVAK) ,pbUF,860) 

K=1 

GMAX=GMAXV( IVAH) 

JVAR=I VAR 

IF  (IVAH-NVAHT)  612*611*611 
L JVAR=6 

IKOWHslHOWT 

> GO  TO  (6125,6125,6121,6122,6122,6125,6122) ,IVAR 
:i  IF  (NSTAT-2)  6125 * 680 , 6i25 

■ 2.  IF  (NSTAl-l)  6125*680,6125 
:5  IF  (LTOLF)  613*614,613 

> IF  (GtfAX)  6132,6132*6135 


I=MMIN(1,LT0PF )»1USV( 1VA«)+1 

WRITE  (NOUT, 1000)  JVCHV(IVAR) , (HDSV(J,I ) , J=1 , 2 ) ,HNUL(K) 

DO  b70  1 YP  = 1 * NYP * 10 

IYPSX=NYPSX+IYP-1 
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IYPEX=MMlN< IYP$X+9,NYP£X) 
n=iypex-iypsx+i 

K=N+1 

NZP=440/K*K 

IF  (LTOLF)  616,618,616 
6lb  IF  (GtfAX)  6l8 « 618  *6165 

6165  WHITE  (NOUT  « 10U3 ) < YOLUV ( 1 ) , 1 = 1 YPSX , I YPEX ) 
WHITE  (NOUT, 1001) 

618  I=u 
U=0 
L=0 

DO  650  NR=NRSTX,NH£NO 
NP(jS=NR 
1=1+1 

CALL  WAKRH  ( NK , Z V ) 

PdUF ( I ) =Z 

IF  (NR-NRST)  6181 <6183,6184 
6l6l  IF  (L2LFF)  6187,6182,6182 

6182  IYPS=IYPS(/(IMR  + 1 ) 

IYPl=IYPEV(NR+i) 

GO  TO  619 

6183  IF  (LZLFF)  6186,6187,6187 

6184  IF  (NR-NREND)  6186,6185,6185 
618b  IYPS=IYPSV(I\,K-1) 

IYPE=IYPEV(NR-1) 

GO  TO  619 

6186  IYPS=KMIN(1YPS,IYHSV(NK-1) ) 
IYPE=PMAX(IYPL,IYHEV(NK-1) ) 

IF  (NR-NREND+1)  6167,619,619 
6167  IYPS=MIWII|\I(IYPS,IYFSV(  NK+1 ) ) 
IYPE=MMAX(IYPL,IYPEV(NH+1) ) 

619  L=L+1 
NN(L)=0 
MM(L)=1 

DU  640  UY=IYPSX,IYPLX 
1 = 1 + 1 
J= J+l 

IF  (JY-IYPS+1)  6195*630 ,620 
6195  HM(L)=MM(L)+1 

GO  TO  640 

620  IF  (JY-IYPE-1)  630,630,621 

621  1=1+1 yPEX-Uy 

j=j+iypex-jy 

GO  TO  641 

63O  CALL  WAKMP(IR0WR,NR*JY,XV,1) 

PdUF(I)=XV< JVAR) 

DBUF ( J ) =XV ( UVArt ) 

NN(L)=NN(L)+1 

640  CONTINUE 

641  IF  (I-NZP)  6415,642,642 
6415  IF  (NR-NREND)  650*642*642 


642  I=I/K 

IF  (LTOPF)  6432,6438,6432 
6432  IF  (IOSV(IVAR))  6434,6436,6434 
6434  CALL  PBFqw ( PlTI D ,NP * U , ObUF ) 
6436  IF  (LTOLF)  6438,646*6438 
6438  IF  ( GMAX ) 646*646*644 


DO  645  JZ=1,I 
IF  (NN(JZ))  645,645*6440 


A-51 


0 

a 

y 

u 

y 

L) 

0 

D 

LI 

0 

0 

0 

0 

li 

* * 

i! 


i - 
• ! 


6440 

L= ( dZ-1 ) *K+ 1 
JS=L+j 

JE=JS+NN ( OZ ) -1 

6441 

GO  TO  (6441,6442.6443,6444.6445.6446.6447,6448,6449,6450) td 

WRITE  ( NOUT , 2031 ) PbUF(L) , (PBUF(d) 
GO  TO  645 

. d=dS , dE ) 

6442 

WRITE  (NOUT, 2032)  PtiUF(L) , (PBUF(d) 
GO  TO  64b 

t d = dS , dE  ) 

6443 

WRITE  (NOUT, 2033)  POUF ( L ) , ( PBUF ( J ) 

. d — d S , dE ) 

GO  TO  645 

6444 

WRITE  (NOUT, 2034)  PBUF(L) , (PBUF(d) 
GO  TO  645 

. d=dS , dE ) 

6445 

WRITE  (NOUT, 2035)  PBUF ( L ),( PBUF ( d ) 
GO  TO  645 

. d — dS , dE  ) 

644fa 

WRITE  (NOUT, 2036)  PBUF ( L ),( PBUF ( d ) 
GO  TO  645 

. d=dS , dE ) 

6447 

WRITE  (NOUT, 2037)  PBUF ( L) .( PBUF ( d ) 
GO  TO  645 

. d=dS . dE ) 

6448 

WRITE  (NOUT, 2036)  PBUF ( L ),( PBUF ( d ) 
GO  TO  645 

. d-dS , dE ) 

6445 

WRITE  (NOUT, 2039)  PBUF ( L ),( PBUF ( d ) 
GO  TO  645 

. d = dS , dE ) 

64bO 

WRITE  (NOUT, 2040)  PBUF ( L) .< PBUF ( d) 

. J=JS, JE) 

64  5 

CONTINUE 

646 

1 = 0 
d = Q 
L=0 

CALL  SFVFL(ZERUVdVAH)  , PBUF  ,880) 

65U 

CONTINUE 

670 

CONTINUE 

630 

CONTINUE 

IF  (LTOHF)  690,700,690 

690 

r 

NRX=NP 

CALL  PBFDW(PLT1D.NRX,1,-1,0) 

C WRITE  COMPLETED  SOLUTION  FOR  NEXT  INTEGRATION  STEP  (WAKSC) 
r 

700 

IF  ( LMLFL ) 702.600,800 

702 

NPTS=NPTSN 
X = XP 

IF  (LSTFL)  705,710,705 

705 

LS>TFL=0 

710 

r 

NRX=1 

CALL  PBFDW(C0M1D«NRX.NC0MT,NSTST) 

C SET  up  for  RETURN  OK  ENU 

r 

600 

IF  (LTRNF)  900,850,900 

850 

N=1 

GO  TO  910 

660 

IF  (KLOCK(l)-MXHRB)  870,660,880 

870 

IF  (dOBE)  875,675.872 

87c 

IF  (dS-dOBE)  875.900.900 

875 

Rt-TURN 

880 

WRITE  (NOUT, 2012) 

900 

N=2 

d5=NG/160+MMIN(l,NlOPP) 

dE=NP/lbO+MMIN(i,NTUHP) 
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WHITE  (NOUT ,20U5)  XP »NPTSN , JS , JE 
910  CALL  CLOCK(l,KLOCH) 

JS=KL0CK(3) 

DO  924  l=lt6.2 

J=b-I 

K= J-2 

JL=KLOH( J)+KL0CK(O) 

IF  (K)  923 » 923 « 921 

921  IF  (JE-59)  923,923,922 

922  Jt= JE-bO 
KLOCK(K)=KLGCK(K)-H 

923  KLOCK(J)=JE 

924  CONTINUE 

GO  TO  ( 860 « 930  ) , N 
930  DO  932  I=l»6«2 

KLOK( I )=KLOcH( 1 ) 

932  CONTINUE 

WRITE  (NOUT.2010)  ( KLOK ( I ) , 1 = 1 ,6 , 2 ) 
NKX  = 1 

CALL  P6F0W<C0MlDtNRX«NC0MT.NSTST) 
CALL  CLbFK(BUFR) 

CALL  CLOFH(BUFS) 

CALL  EXIT 
END 
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♦♦wakpc  - stratified  submarine  wake,  pressure  calculation 
SUBROUTINE  wakpc 
c 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  CONTROLS  THE  SOLUTION 
C FOR  THE  PRESSURE  BY  THE  POISSON  EQUATION 
C 

DIMENSION  TOTAL  (7, 40, 2)  ,PN(  5 ) , V ALUE  ( 40 , 2 ) 

♦COPY  (CMwAK) 

F.UU1  VALENCE  (TOTAL  (1,1,1)  ,KOWG(l,l)  ) , (MM,EPSXV(1)  ) 

LUOI VALENCE  ( ERROR, bM( 12) ) , ( GM ( 10 ) , DUI ) , ( GM ( 9 ) , OLI ) 
EQUIVALENCE  (VALUE (1,1) ,ROWG(l,36) ) 


C 

DATA  pN/iO,  ,30*  ,6’J,  *100  , , 1 5 0 « / 

C 

1000  F0RMAT(//6X,6HSTATUS,7X,SHXX  PO INTS , bX , 9HMAX  VALUE , bX , 

1 lOhMAX  cHANGE»8X, ‘iHUXX,HX,5HERR0R/8X,2HPR,  I14«5X'/4E15.5) 
2010  FORMAT ( //41H  PROGRAM  EXIT  LURING  PRESSURE  COMPUTATION) 

C 

C CHECK  RUN  POSITION 
C 

I OLAY =8 

IF  (LSTFL)  SO, 105*100 
50  LSTFLrLPRFL 

IF  (LSTFL)  310,65*8b 
85  MPTS=MM 

DXXrEPSXV(2) 

GO  TO  109 

100  IF  (LFRFL)  310, 10b, 105 
10b  MPTS=0 

NREND=NRST+JEND 

Y=YOLDV ( NYPE  + 1 ) -YOLUV ( NYPS-1 ) 

Z = ZOLDV.(  NKEND  ) -ZOLDV  ( NRST-1 ) 

DXX=PNORM+Y+Z/PN(l) 

109  LPRFL=-1 
MOOu=0 
IBOT=l 
IT0P=MXRY 
LYFAFsO 

CALL  SFVMV(ZEROV,XPHV,NWVEC) 

CALL  sfvfl(o.o,zerov,nwvec) 

IF  (LFOLF)  1090,110*1090 
1090  NYPSX=NYPS+lYlFF 

00  1094  JY=NYPSx,N'/PE 
IF  (YOUT-YOLOV(JY) ) 110,1096,1092 
1092  IF  ( YOUT-YOLOV(UY+l) ) 1096*1094*1094 
1094  CONTINUE 
1096  LYFAFrJY 

RATY=( YOUT- YOLUV ( 0 Y > ) / ( YOLOV ( JY+1 ) -YOLDV ( Jy ) ) 

CALL  SFVFHO.O, VALUE, 80) 


|l 

L 


c 

C INITIALIZE  FOR  NEXT  INTEGRATION 
C 


CALL  OATSW ( 0 , ISW ) 

GO  TO  (111, 114), ISW 
WRITE  (NOUT ,2010 ) 
MM=MPTS 
EPSXV ( 2 ) =OXX 
LPRFL=LSTFL 


STEP 


I 

I 

] 
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X=XP 

NHTSsNPTSN 
GO  TO  300 
U*»  0X1=2, O/OXX 

FMAXV ( 7 ) =0 , 0 
TMAXV ( 7) =0 , 0 
OUIsO.O 
DLI*0,0 
MPTS=MPTS+1 
CALL  WAKPZ 
C 

C DETERMINE  NEW  OXX  FOR  PRESSURE  SOLUTION 
C 

QXXOsQXX 

IF  (MPTS-NSTPR)  127*  170*170 
127  OXX=DXX*PN(MPTS}/PN(MPTS+l) 

C 

C SUCCESSFUL  INTEGRATION  SjEP 
C 

170  EKKOR=SQRT(OUI/OLl )/OXXO 

WRITE  (NOUT , 1000 ) MPTS , FMAXV ( 7 ) . TMAXV ( 7 ) ,DXXO , ERROR 
IF  (MPTS-1)  172,172*171 

171  IF  ( TMAXV ( 7 ) -TEM ) 172,172,180 

172  TEMsTMAXV ( 7 ) 

IF  (MPTS-NSTPR)  175*180,180 
175  IF  (TMAXV(7)/FMaXV(7)-PCRIT)  180. 180*110 

180  IF  (LFOLF)  210,29U,210 

210  NRSTXsNRST+LZLFF 

nreno=nrst+jeno-i 

00  230  NR=NRSTX,NKENO 
IF  (ZOUT-ZOLQV(NR) ) 285*240,220 
2:0  IF  <ZOUT-ZOLOV<NR+l) ) 240.230,230 
230  CONTINUE 

240  RATZ=(ZOUT-ZOLOV(NR) ) / < ZOLDV < NR+1 ) -ZOLDV (NR ) ) 

NPOSsNR 

CALL  WAKRR(NK.ZV) 

CALL  WAKRR(NR+1,ZPV) 

IYPSXsIYPS+LYLFF 

00  260  UY=IYPSX,IYPt 

CALL  WAKMP(IROwT,NR, JY,XV,1J 

CALL  WAKMP(IROWT,NR+l,JY,XPYV,l) 

T0TAL(7,JY,1)=XV(6)^KATZ*(XPYV(6)-XV(6) ) 

260  CONTINUE 

265  IF  (LYFAF)  290,290,270 

270  DO  280  NR=NRSTX,NKENU 

TOTAL ( 7 vNR*2 ) =V ALOE (NR, 1 ) +K ATY *( VALUE (NR,2) -VALUE (NR, 1) ) 
280  CONTINUE 
290  GMAXV ( 7 ) =FMAXV ( 7 ) 

300  CALL  SFVMV(XPPV,ZEROV,NWVEC) 

310  RETURN 
END 
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WARP  I . S ( 01 05  ) 

**wakpi  - stratified  submarine  wake 
subkoutine  wakpi 


profile  INITIALIZATION 


c 

C This  SUBROUTINE  IN  THE  WAKE  PROGRAM  INITIALIZES 
C THE  WORKING  FILE  WITH 

c 

C 1)  KNOWN  JET  PLOT  FILE  RESULTS  FOR  AXISYMMETRlC  FLOW 
C 

C 2)  KNOWN  STRATIFIED  WaKE  PLOT  FILE  RESULTS  FOR  SPECIAL  RESTART 
C 


0 

fl 

0 

0 


n 

o 


0 


C 3)  GIVEN  PRESSURE  DISTRIBUTIONS  FOR  POISSONS  EQUATION 

c 

DIMENSION  FILlDt  3> * RBUF ( b , 4 0 ) tRH( lb) , JT0P(2) 

♦COPY  (CMWAK) 

EQUIVALENCE  (KoUF(ltl) ,ROWG(l,l) > ♦ ( A V ( 1 ) . F IL I D < 1 ) ) 
EQUIVALENCE  ( J10P<2) , ZNEW V ( 2 ) ) , ( RR ( l ) , GM ( 1 ) ) 

C 

1000  FORMAT ( //24H  IMPROPER  INI  1 1ALIZATI0N.4I5) 

C 

I0LAY=4 

CALL  SFVFL(O.OtGMAXVtNVART) 

MOOD=-l 

IF  (LyFAF)  210*200,50 
C 

C FIND  APPROPRIATE  JET  PROFILES  TO  INITIALIZE  NSTAT  = 1 
C 

50  NR=2 

NVARI=3 
10U  NP=NK-1 

CALL  PBFDR(FILIDtNPtLZFAF,ZNEWV) 

I T OPsjTOP ( 1 ) 

NR=NP+ITOP*LYFAF 

CALL  PBFDRtFlLID.NRfl.XP) 

IF  (XP)  103  1 102  1 102 

102  IF  (XP-X)  100*100*103 

103  NR=1 
X=ZNEWV(1) 

CALL  SFVFL(0.O«XV*NVAR) 

XV ( 6 ) =SC ALE 

IBOT=MMIN(MxRY/(LTLFF+2) «MXRZ/(LZLFF+2) ) 
NSK=ITOP/IBOT+l 
JO  = 0 
NY  = 0 

DO  120  J=1 * I TOP , NSK 
NY=NY+1 

NP=NP+ ( J- JO-1 ) *LYF AF 
JO= J 

CALL  PBF0R(F ILID*NP«LYFAF*RR  ) 

RBUF ( 1 . NY ) =RR ( 5 ) +RR < 6 ) + RR ( 7 ) 

RBUF(2«NYi=0.0 
RBUF ( 3 » NY ) =RR ( 3 ) 

RBUF ( 4 « NY ) =RR ( 1 ) 

00  115  I=1«NVARI 
TEMsABS ( RBUF ( I « NY ) ) 

IF  ( TEM-GMAXV( i ) ) 115*115.114 

114  GMAXV(I)=TEM 

115  CONTINUE 

IF  (J-ITOP)  lib,  120*120 
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6 IF  (J+NSK-ITOP)  120*120,118 
8 J=IT0P-NSK 
0 CONTINUE 

CALL  WAKSEiEPSS'EPSSV) 

2 CALL  SFVFL(0.0«RBUF(1«NY) tNVARl) 

00  125  Isl'NVARl 

IF  <ABS(RBUF(I*NY-1) )-£PSSV(I)  ) 125,125,126 

5 CONTINUE 
NYsNY-1 
60  TO  122 

6 IF  (L2LFF)  1262,1264,1264 
62  NRST=2 

NRSTXzl 
JEND=NY-2 
GO  TO  1266 

64  NKST=(MXRZ+1 ) /2-NY+2 
NHSTX=MXRZ/2 
JLN0=NY+NY-3 

66  IF  (LyLFF)  1268,1270,1270 
66  NYPS=2 
NYPSXzl 
NTPE=NY-1 
GO  TO  127 

70  NYPS=(MXRY+l)/2-NY+2 
NYPSXsMXRY/2 
NYPE=NYPS+NY+NY-4 

7 IF  (JEND)  129*129,1275 

7 5 IF  (NRST-2)  129,128*128 
» IF  (NYPS-NYPE)  1285*129,129 
ib  IF  (NYPS-2)  129,130,130 
? WRITE  ( NOUT « 1000 ) NYPS,NYPE ,NRST, JEND 
CAL.  EXIT 

CONSTRUCT  AXISYMMETR1C  WOHMNG  FILE 

) DO  180  JZ=1,NY 
NK=NRSTX+jZ-1 
NPOS=NK 
Z=RBUF(4, JZ) 

IF  (JZ-NY)  132*131,131 
L IYPSzo 
IYPEsO 
GO  TO  170 
: J=JZ 

iyps=ntpsx-lylff 

00  160  JY=1,NY 
J=J-1 

NRY=NYPSX+JY-1 
Y=RBUF(4, JY) 

IF  (JZ-1)  133,133,135 
* YOLOV ( NRY ) =Y 

IF  (LYLFF)  135,134,134 
I Y=NYPSX« JY+1 
YOLOV( IY>=-Y 
R=SQR7( Y*Y+Z*Z) 

J=J+1 

IF  (JY-NY)  141,165,165 
IF  (R-RBUF (4, J+l ) ) 142,142,140 
RATR=(K-RBUF(4, J) )/<KBUF< 4, J+l )»RBUf(4, J) ) 


«■ 
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c 

C F 
C 

<200 


CAKT 


CALL  SFVMV(RB'JF(1«J)  tXZM'SNVARI) 

CALL  SFVMVIRBUF  <1« J+l) iXZPV,NVARI ) 

I VAR  = Q 

00  146  1 = 1 «f\IVARI 

XV  (I  )=XZMV(I  )+RATR*».'<ZPV(I)-XZM\/(  I)  ) 

IF  (ABS(XVd)  )-EP5»SV(I)  ) I45*14b»146 

IVAR=IVAR+1 

CONTINUE 

IF  (IVAR.NVARI)  lb0«165t)bb 
CALL  WAKMP(IR0wR,NR.NRY,XV«2) 

IF  (LYLFF)  I60«lbl,lb). 

IF  (JY-1)  160»lfa0«lS2 
IY=NYPSX-JY+1 

CALL  WAKMP( IR0wR«NR»lY»XVf2) 

1 YPS=I YPS-1 
CONTINUE 
IYPE=NRY-l 
CALL  wakwk;nr*zv) 

IF  (L2LFF)  i80  » 171 » 171 
IF  (JZ-l)  180tlB0«l72 
NR=NRSTX-JZ+1 
NPOS=NR 

z=-z 

CALL  UAKMK(NR«^V) 

continue 
GO  TO  2b0 

1ND  APPROPRIATE  SPECIAL  RESTART  PROFILES 


CALL  WAKSR(FILIQ) 

GO  TO  2b0 

CALL  WAKSK(PLTiO) 

NP=1 

NG=1 

RETURN 

ENU 
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WAKPM, S ( 0105 ) 

♦♦WAKPM  - STRATIFIED  SUBMARINE  BAKE*  PRESSURE  MANIPULATIONS 
SUBROUTINE  WAKPM(LFL) 


1105 


1205 

1208 

1209 

121 

1215 

C 

122 

123 

C 

130 


THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  PERFORMS  VALUE  Ma( 
ON  THE  UPSWEEP  AND  OCWNSwEEP  OF  THE  PRESSURE  IMPLICIT 

DIMENSION  VALUE<40,2) 

♦COPY  (CMWAK) 

EQUIVALENCE  (GM(10) tOUI ) « (GM(9) , DLI ) 

EQUIVALENCE  (VALUL(1,1) , K0W6 (1*36) ) 

SPRZ( ARGM* ARG, ARGP ) =SZM*ARGM+SZ*ARG+SZP*ARGP 
SPRY( ARGM* ARG»ARGH )=SYM*ARGM+SY*ARG+SYP*ARGP 

GO  TO  (100,200) ,LPL 

matrix  upsweep  CALCULATIONS 

) OZMsZ-ZM 
DZPsZP-Z 

dzt=dzm+dzp 

SZMs2,0/DZM/DZT 
SZP=2,0/DZP/0ZT 
SZr-SZM-SZP 

dym=y-ym 

DYP=YP-Y 
DTTsDYM+DYP 
SYM=2 ,0/DYM/DYT 
SYP=2.0/DYP/0YT 
SY=-SYM-SYP 
GO  TO  (110,120) ,LZFAF 

I XMAT(1)=-SYM 
YMAT(1)=DXI-SY 
ZMAT(1)=-SYP 

DVEC(1)=DXI*XV(6)+SPKZ(XMYV(6),XV(6) . XPYV ( 6 ) ) -TV ( 5 ) 
IF  (JY-2)  1105,112,111 
15  TLM=FLOAT(IRFV<6,3)  ) 

GO  TO  122 

IF  (JY-IYPSN)  123*123,130 
IF  (LYLFF)  130,123,123 

XMAT(1)=-SZM 
YMAT(l)=OXI-SZ 
ZMAT ( 1 ) =-SZP 

DVEC(1)=DXI*XV(6)+SPKY(XZMV(6),XV(6) ,XZPV(6))-TV(5) 
IP  (JZ-1)  1205,1215*1206 

5 TEM=FlOAT (IRFV(b,H) ) 

GO  TO  122 

6 IF  (JY-IYPSM)  1209,130*121 
9 IF  (JY-1)  130*130*123 

IF  (JY-IYPEM)  130*130,123 
5 IF  (LZLFF)  130,123,123 

ZMAT(l)=ZMAT(l)  + (2.0nEM-1.0)^XMAT(l) 

XMAT ( 1 ) =0 • 0 

TEM=l,0/(YMAT(l)-XMAT(i)*AV(l) ) 


value  manipulations 

IMPLICIT  SOLUTION 
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AVtl»=TEM*ZMAH  1) 

AV(b)=TEM*lOVEC(l > -XMAT < 1 ) *AV (6 ) ) 
CALL  WAKMP(IR0*A,NR*JY,AV,2) 

RlTURN 

c 

c matrix  oownsweep  calculations  <wakpy> 
c 

200  CALL  WAKMP< 1R0WA,NR*JY,XV,1) 
XV(6)=XV(6)-XV(1)*AV<6> 

A V ( 6 ) zXV ( 6 ) 

CALL  WAKMPtIRO*T,NR,JY,TV,l) 
TLM=ABS(AV(6)  ) 

IF  ( TEM-FMAXV ( 7 ) ) 240,230,220 
220  FMAXV(7)=TEM 

YMAXV(7)=YOLOVl JY ) 
ZMAXV(7)=20LDV(NR) 

230  TLM=ABS(AV(6)-TV(6)  > 

Y=YULDV( JY ) 

Z=ZULDV(NR) 

DY=YOLOVldY+l)-Y 
IF  (JY-1)  233,233,242 
232  UY=DY+Y-YOLDV( JY-1 ) 

234  OZ=ZOlOV(NR+l)-Z 

IF  ( NR-1 ) 235,235,244 

234  OZ=DZ+Z-ZOLUV(NR-l ) 

235  DUI=bUl+TLM+TEM*OY*UZ 
0LIrQLl  + TV(5)*l  V(b)*UY*QZ 
IF  ( TEM-TMAXV ( 7 ) ) 250,250,240 

240  TMAXV ( 7 ) =TEM 

250  TV  ( 6 ) =AV ( 6 ) 

CALL  WAKMP( IROwT.NR, JY,TV,2) 

IF  (LYFAF)  260,260,252 
252  IF  (JY-LYFAF)  260*254,256 
254  VALUE(NK,1)=TV(6) 

return 

256  IF  (JY-LYFAF+1)  260*250,260 
256  VALUE(NK,2)=TV(6) 

260  RETURN 
ENb 

CAKT  ID  0105  OB  ADDR  4240  DB  CNT  0QC2 
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WAKPP.S< 0105 ) 

**WAKPP  - stratified  submarine  wake,  printer  plotter  routine 
subroutine  wakpp 
c 

c THIS  SUBROUTINE  IN  THE  WA«E  PROGRAM  SCANS  THE  WORKING  FILE 
c aNO  CONSTRUCTS  a point  plot  fur  the  printer 
c 

DIMENSION  11(21) ,UVtC(101) 

DIMENSION  XVEC(40>  ,UVCHV(7)  * JSCAL(lOl) 

♦COPY  ( CMWAK ) 

EQUIVALENCE  ( JVEC ( 2 ) * YNEW V ( 1 ) ) « ( XVEC ( 1 ) , ZNEWV ( 1 ) ) 

C 

DATA  IH/1H-/«IV/1HI/«IO/1HO/ 

DATA  II/lHM«lHH'lh  * 1H8  * 1H  ,1H6.1H  ,1H4,1H  ,1H2,1H  ,lHl,iH  , 

1 1H3,1H  ,1H5,1H  « 1H7 « 1H  ,lHp,lHP/ 

DATA  jVCHV/2HQQ,2HnO,2HU  ,2HV  , 2HW  ,2HSL,2HP  / 

C 

1005  FORMAT ( 1H1 ) 

1010  FORMAT ( 18H  PRINTER  PLOT  FOR  *a2.5X,11HX  POINTS  = ,I3,5*,4HX  s , 
1 E11.5.5X.16HABS  MAX  VALUE  = ,E11.5,5X,13HEDGE  VALUE  = ,E12,5) 
1020  FORMAT(20X,101A1) 

1040  FORMAT (12H0PERCENT  MAX ,3X , 6H90/100 ,5X, 5H70/80 ,5X , 5H5G/60 , 5X , 

1 5h30/40,5X,5H10/20«5X,4H1/-1 .5X*7H-10/-20,3X»7H-30/-40,3X» 

2 7H-50/-60,3X,7h-70/-80,2X*8h-90/-100/12H  NOTATION  *bX,lHP, 

5 10X,1H7,9X,1H5,9X,1H3,9X.1H1,8X.1H0*10X*1H2*9X.1H4,9X*1H6*9X, 

4 lH8t8X.lHM/) 

lObO  FORMAT(E18.5«1X,103A1) 
lObO  FORMAT(l9X.103Al) 

1070  FORMAT ( 6X  *6E20  • 5 ) 

1060  FORMAT(5x.1HZ.13X*103A1) 

1090  FORMAT(/70X,1HY) 

C 

C INITIALIZE  PLOT  VALUES 
C 

IOLAYz9 

IF  (LPRFL)  I0«350*3b0 
10  IF  (LTOLF)  15*20.15 
15  IF  (LOUT)  20*50*30 
20  CAUL  DATSW(7«ISW) 

GO  TO  (bo « 350 ) « ISW 
3o  CALL  DATSW(3«ISW) 

GO  TO  (bO  *40 ) « ISW 
40  LTOLF=0 

50  MOOD=l 

NY=101 
NZ=51 

NYPSX=NYPS-1 

nypex=nype+i 

nrstx=nrst-i 

nrend=nrst+jenu 

c 

C SET  UP  AXES  AND  NOTAIION 
C 

YMlNrYOLOV(NYPSX) 

YMAXsYOLOV(NYPEX) 

ZMINzZOLDV(NRSTx) 

ZMAXzZOLDV(NRENO) 

DYs(YMAX-YMIN) /FLOAT (NY-1) 

DZ= (ZMAX-ZMIN)/FLOAT(NZ-l) 
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0 

II 


U 

(J 

IJ 


00  80  1=1,6 

YMAT( I)=YMIN+20.G*FLOAT( I-1)*DY 
2MAT( I )=ZMIN+10.0*FLOAT(I-1)*DZ 
80  CONTINUE 
K=1 9 

00  90  1=1, NY 
K=K  + 1 

IF  (K-20)  84,82,82 
82  K = 0 

JSCAL(I)=IV 
GO  TO  90 

84  JSC AL ( I ) =IH 

90  CONTINUE 

C 

C LOOP  THROUGH  FILE  FOR  EACH  VARIABLE 
C 

IKOWR=2 

IF  (LBUKF)  92*94,92 
92  IK0WR=4 

94  WRITE  ( NQUT « 1005 ) 

DO  300  IVAR=1,NVAKT 
FHAX=FMAXV(IVAR) 

FMIN=l.0C-04*FMAX 
IF  ( FMAX ) 300,300,100 

100  JVAR=IVAr 

IF  (IVAK-NVART)  .02*101, 101 

101  JVAR=6 
IROWR=IROWT 

102  WRITE  (NOUT,lOlO)  JVCHV(IVAR) » NPTSN , XP , FMAX , ZEROV(IVaR) 
WRITE  ( NOUT  1 1040 ) 

WRITE  (NOUT, 1020)  JSCAL 

KS=9 

IS=7 

NS=NZ/2+l 

nr=nrenu-i 

npos=nr 

CALL  WA)RR(NK+1 ,2PV) 

npos=nr+i 

CALL  WAKRR ( NR , ZMV ) 

C 

C LOCATE  EVERY  OESIREO  2 VALUE 
C 

DO  250  12=1, NZ 
JZ=NZ- IZ+1 

Z=ZMIN+FL0AT(J2-1 )*02 

105  IF  (2-2P)  107,110*106 

106  Z=ZP 

GO  TO  110 

107  IF  (Z-ZM)  108*110*110 

108  IF  (NR-NRSTX)  1095,1095,109 

109  NR=NR-1 

CALL  WAKMR (1,3) 

CALL  SFVMV(Z«,ZP,NWWZf 

NPOS=NR+l 

CALL  WAKRR ( NR , ZMV ) 

GO  TO  105 
1095  Z=ZM 
C 

C interpolate  FOR  PaRaCETER  VALUES 
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C 

lie  RATZ=(Z-ZM)/(ZP-Zfl) 

00  115  JY=NYPSX,NYPEX 

CALL  WAKMP(IR0*R,NR,JY,XMYV,1) 

CALL  WAKMP(IR0*R,NR+2,JY,XPYV,1) 

XV£C(  JY)=XMYV(JVAH)+RATZ*(XpYV(  JVAR)-XMY''(dVAR)  ) 
lib  CONTINUE 
C 

C EXPAND  TO  FILL  PLOT  ARRAY 
C 

K=NYPSX 

DO  160  JY=liNY 
Y=YMIN+FLOAT( JY-1)*UY 

116  IF  ( Y-YOL0V ( K ) ) 119*120,117 

117  IF  (Y-YOLDV(K+l) ) 120 ,120,118 
lib  K=K+1 

IF  (K-NYPEX)  116*119*119 

119  TEM=XVEC(K) 

GO  TO  121 

120  RATY=( Y-yOLDV(K) )/( rOLOV(K+l)-yOLOV(K) ) 
TEM=XVEC(K)+RATY*(XVEC(K+1)-XVEC(K) ) 

c 

C CHECK  AND  INSERT  POINT  VALUES 
C- 

121  I = IFlX(10.0*TEM/FriAX) 

JVEC(JY)=II(I+11) 

IF  (I)  160*122*160 

122  I=IFIX(100.0*TEM/FMAX) 

IF  (I)  160*124*160 

124  IF  (ABS(TEM)-FMIN)  160*126*126 
12&  JVEC ( JY ) =10 
160  CONTINUE 
C 

C print  TO  printer 

C 

ks=ks+i 

IF  (KS-10)  204,202,202 
202  KS- 0 

IS=IS-1 

WRITE  (NOUT , 1050 ) Z^AT < IS ) * IH , JVEC , IH 
GO  TO  2b0 

204  IF  (JZ-NS)  208*206,208 
206  WRITE  (NOUT, 1060)  IV,JVEC*IV 

GO  TO  2b0 

206  WRITE  ( NOUT ,1060)  IV, JVEC, IV 

250  CONTINUE 

WRITE  (NOUT, 1020)  JbCAL 

WHITE  (NOUT, 1070)  ( YMaT ( I ) * 1 = 1 , 6 ) 

WRITE  (NOUT, 1090) 

300  CONTINUE 
350  RETURN 
END 

ca«t  id  0105  00  aoor  4eco  ob  lnt  0152 
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WAKPY.S(0105) 

**W AKPY  - STRATIFIED  SUBMARINE  WAKE,  PRESSURE  SI  EP  IN  Y 
subroutine  WAKPY 


THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  MAKES  A STEP  OF  DELTA  XX/2 
WITH  i DERIVATIVES  EVALUATED  AT  THE  PRESENT  X*  POINT 


*CUPY  (CMwAK) 

C 

C INITIALIZE  COMPUTATION 


I YPSNzIYPS 
IYPtN=IYPt 


INITIALIZE  FOR  THE  Y DIRECTION  AND  SWEEP 


lObl 


I0b2 

1063 


Y=YOLDV< IYPSN-1 ) 

YPrYOLOV(IYpSN) 

CALL  WAKH?( IR0WR,NR+1 , I YPSN-1,XPYV*1 ) 
CALL  WAKMP(lROwR,NR, 1 YPSN- 1 , X V , 1 ) 

CALL  WAKMP(IROwT,NR,IYPSN-l,TV,l> 

IF  (JZ)  1081,1081,1042 
CALL  wAKRF(XPYV,XI“YV,4) 

GU  TO  1085 

CALL  WAKMP(IR0WK,NR-1, IYPSN-1, XMYV»1) 
CALL  WAKMP( IR0WR,NR, IYPSN,XZPV,1) 

JY=1 YPSN-1 

IF  (LyLFF)  109,110,110 
CALL  WAKRF(XZPV,XZMV,3) 
ym=y+y-yp 
GO  TO  120 


UPWARD  PASS 


JYrJY+1 

CALL  SFVMV(XV,XZMV,NMOVE) 

CALL  WAK«P(lROwT,NR, JY,TV,1) 

YM=Y 
Y = YP 

YP=YOLDV( JY+1) 

CALL  WAKMP(IROwR,NR'fl,JY,XPYV,l) 
IF  (JZ)  116,116,118 
CALL  WAKRF  (XPYv  , X P*  Y V , 4 ) 

GU  TO  119 

CALL  WAKMP(IROwrt,NR-l , JY,XMYV, 1 ) 
CALL  WAKMP(lROwR,NR,JY+l,XZPV,l) 


calculate  Matrix  coefficients  and  invert  for  gamma  and  av 


CALL  VAKPM(I) 

IF  (JY-IYPEN)  110*126, 125 


UPPER  BOUNDARY  CONDI  I ION 


TV(6)=AV(6) 

CALL  WAKMP(IR0mT,NR,JY,TV,2) 


INITIALIZE  FOR  uOwNWARU  PASS 
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150 


IYPENtJY 
I YP£.y.=  7YP£N  »1 
i YKiaXsI  YPSN+LYLFF 
IF  (IYPSX-IYPCN)  150  « 156 1 156 
DO  151  IY  = IYPSX,IYP£.X 

jy=iypex+iypsx-iy 

CALL  WAKPM( 2 ) 

CONTINUE 


151 
C 

C ROw  SWEEP  COMPLETED 
C 

15b  IYPS=IYPSN 
IYPE=IYP£N 
CAUL  WAKWR(NRt^V) 
RETURN 
END 

CART  10  0105  OB  AODR  3230 
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WAKPi.S(Q105> 

**WAKPZ  - STRATIFIED  SUBMARINE  WAKE,  pressure  step  in  Z 
subroutine  wakpz 
c 

c this  SUBROUTINE  in  the  wake  program  MAKES  a STEP  of  delta  XX/ 2 
C WITH  Y DERIVATIVES  EVALUATED  AT  THE  PRESENT  XX  POINT 
C 

♦COPY  (CMWAK) 

C 

C INITIALIZE  COMPUTATION 

c 

LZFAF=2 

IROwRsl 

IKOWAri 

nrstn=nrst 

jenun=jeno 

NR=NRsTN-1 

npos=nr 

CALL  WAKrR(nK.ZV) 

CALL  WAKRR(NR+i,2PV) 

JZ  = 0 

IF  (LZLFF)  lOl  * 10 15 1 10 15 

101  Z M=Z+Z-ZP 

I YPEMrl YPEP 
IYPSM=IYPSP 
GO  TO  llo 
1015  NRrNRSTN 

102  JZ=JZ+1 
NPUS=NR 

CALL  SFVMV(ZtZM,NWW2F) 

CALL  WAKMR(2tl) 

CALL  SFVMV(ZP»Z,NWWZF) 

CALL  WAKMR (5*2) 

CALL  WAKRR ( NR+1 1 ZPV ) 

C 

C initialize  FOR  THE  y direction  and  sweep 

C 

HO  Y = YOLDV(  IYPS-1) 

YP=YOLUV( IYPS) 

CALL  WAKMP(IROWR,NR+ltIYPS-ltXPYVil) 

CALL  WAKMP(lRO*R, NR, IYPS-1, XV, 1) 

CALL  WAKMP ( IRQWT , NR, IYPS-1, TV, 1) 

IF  (JZ)  111. 111. 1HS 

111  CALL  WAKRF(XPYV,XRYV.4) 

GO  TO  112 

1115  CALL  WAKMP(IROwK,NR-l, IYPS-1, XMYV.l) 

112  IYPSXrlYPS+LYLFF 

IF  (JZ)  1126.1126.1125 
112b  CALL  WAKMP(IROwA.NR-l.IYPS-l.AV.l) 

1126  CALL  WAKMP(IR0«R,NR.XYPS»XZPV,1) 

DO  145  Uy=IYPSx,IYPL 
IF  (JY-IYPS)  114.1142,1142 
114  CALL  WAKRF(XZPV,XZMV,3) 

YM=Y+Y- YP 
GO  TO  120 

11^2  IF  (JZ)  1145,1145,1144 

1144  CALL  WAKMP( IROwA.NR-1, JY.AV.l  ) 

1145  CALL  SFVMV(XV,XZMV,NMOVE) 

CALL  WAKMP( IROwT.NR, JY.TV.l) 
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ii4e 


1149 

L 


YM=Y 
Y = YP 

YP=Y0l0V< JY+1) 

CALL  WAKMP(lR0wR,NRU,JY,XPYV,l> 
IF  (J2)  1148,1148,1149 
LALu  WAKRF(XPYV,XKYV,4) 

GO  TO  119 

CALL  WAKMP< IRO»R,NR-l,JY,XMYV,l) 
(.ALL  w AKMP  ( lROwR,NR,JY+l,XZPV,l) 


calculate  matrix  cou-ficilnts  ano  invert  for  gamma  ano  a v 


CALL  WAKPM(i) 

CONTINUE 

IF  ( JZ-JEnDN+1 ) 146*147,147 
CALL  wAKwR(NR,Zv) 

NK=NR+1 

IF  (J2-UENDN)  102*150,150 


initialize  fur  downward  pass 


1502 

1503 
150s 


NK=NR-1 

JENDN=JZ 

izend=jendn-lzlff-i 

LZFAF=1 
IR0WRs3 
Irt0wA.=  4 

DO  154  1 2 = 1 « IZEND 
iyend=iype-iyps-lylff+i 

DO  1508  I Y = 1 * I Y END 
JY=1YFE-IY+1 

IF  (JY-IYPEM)  1502,1502,1508 
IF  (JY-IYPSM)  1505 | 1506 ,1506 
IF  (IYPSM-2)  1506,1506,1506 
CALL  WAKMP( 1ROwR,NR, JY,AV,1  ) 
CALL  UAKMP(lROrtR,NR-l,JY,XV,l) 
XV(6)=XV(6)-XV(1)*AV(6) 

CALL  UAKMP( IR0WR,NR-1, JY,XV  ,2) 
CONTINUE 


1508  CONTINUE 
C 

c row  sweep  completed 
c 

CALL  WAKPY 
NR=NR-1 

npos=nR 

JZ=JZ-1 

CALL  SFVMV (2,ZP,NWW2F ) 
CALL  WAKNR<2,3) 

CALL  SFVMV(ZM,Z,NWWZF) 
CALL  W AKMR (1,2) 

IF  (IZ-I2END)  153*152,152 

152  IF  (LZLFF)  154,155,153 

153  CALL  WAKRR(NR-1,ZMV) 

] 54  CONTINUE 

IF  (JZ)  156,156,158 
156  ZM=Z+Z-Zp 

iypem=iypep 

IYPSM=IYPSP 
158  CALL  WAKPY 
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WaKRF , S ( 01 05 ) 

**UAKRF  - STRATIFIED  SUBMARINE  WAKE*  REFLECT  A POINT 
SUBROUTINE  WAKRF ( I EMAtTEMBtLFL ) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  PROVIDES  REFLECTION 
C 

DIMENSION  T£MA(b) iTLMB(b) 

♦ COPY  ( CMWAK ) 

C 

CALL  SF\/MV(T£MAtTtMBtNWVEC) 

DU  100  I=ltNWVEC 
IF  C I RF V C I • LFL ) ) I00t50,100 
50  TEMB ( I )=-TEMb( 1 ) 

GO  TO  (6C»60»100tl0U)  ,LFL 
60  TEMb(I)=TEMB(I)+2.0*ZEROV(l) 

100  CONTINUE 
RETURN 
ENu 

Cart  ID  0105  OB  aUDR  4A20  OB  CNT  0022 
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WA*RK.S(0105> 

**WaKRR  - STRATIFIED  f-oMARINE  WAKE*  READ  A Z BUFFER  ROW 
SUBROUTINE  WAKH  ‘l\RX,ZP0S) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  READS  A i ROW  AND 
C FILLS  THE  ROW  BUFFER  AT  THE  NRX  POSITION  FOR  ALL  Y 
C 

DIMENSION  ZP0S(2) «RLCB(24) 

♦COPY  (CMWAK) 

EUUIVaLEN^E  (RECB(l) »DVEC(1) ) 

C 

DATA  JERRX/2HRR/ 

C 

C FILL  ZPOS  WITH  IYPS  ANU  IYPt  VALUES 
C 

IF  (NRX)  300*300*10 
10  IF  (NRX-MXRZ)  20,20*300 
20  ZA=ZOl_DV(NRX) 

IVECA-IYPSV(NRX) 

IVECBrlYPEV(NRX) 

CALL  SFVMVlZA«ZPOS«NWWZF) 

C 

C LOCATE  NRX  POSITION  IN  ROW  BUFFER  ARRAY 
C 

J=NKX-NP0S+2 
IF  (J)  300*300,30 
30  IF  (J-3)  <4  0*40*300 

40  IF  (IVECU  45*45*50 

45  I YPSXsl 

GO  TO  140 

50  IYPSX=IVECA+LYLFF 

IF  (IYPSX-I)  90,90,60 
60  I YPEXsI YpSX-1 

C 

C FILL  BUFFER  ROw  WITH  LLAOING  ZEROES 
C 

DO  60  1=1 , IYPEX 

CALL  SFVFL(0.0»ROWB(1*I*U) *NWVEC) 

CALL  SFVMV(ZER0V»H0WB(7*I* J) »NWVEC) 

CALL  SFVMV(ZEK0V»R0WB(13*I»J) »NWVEC) 

80  CONTINUE 

90  IYPEX=IVECB 

C 

C FILL  bUFFER  RO*  WITH  NRX  INFORMATION 
C 

NRXX=(NRX-1)*MXRY+IYPSX 

DO  130  IrlYPSX, IYPEX 

CALL  PBFDR(SLNID*NRXX,NWR*KECB> 

CALL  SFVMV ( RECB ( 1 ) ,KOWB(l«I*J) ,NWVEC) 

IF  (MOOU)  105*110  *113 

105  CALL  SFVMV(R£CB(7) »RUw8(7*I*U) ,NMOVE) 

GO  TO  130 

110  CALL  SFVMV(RECb(13) »ROWB(7tI»J) »NMOVE) 

GO  TO  130 

115  CALL  SFVMV(RECa(7) ,R0WB(7*I,U) ,NWVEC) 

CALL  SF VMV ( RECB ( 19 ) »R0WB( 13,1 ,J) ,NWVEC) 

130  CONTINUE 

IYPSX=I YPEX+1 

IF  (IYPSX-MXRY)  140*140*160 
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FILL  bUFFER  ROW  WITH  TRAILING  ZEROES 

140  00  150  I=IYPSX,MXHY 

CALL  SFVFL(U*0«ROWQ(1,1,J) «NWVEC) 
CALL  SFVMV(ZEK0V.R0Wd(7.I»W) ,NWVEC) 
CALL  SFVMV(ZEK0VfR0WB(13»I«J) .NWVEC) 
150  CONTINUE 
160  RETURN 
300  JERR= JERRX 
RETURN 
ENO 

CART  iD  0105  DB  AOOR  4A5O  OB  CNT  008C 
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WAKSC.S(Q105> 

♦♦WAKSC  - STRATIFIED  SUBMARINE  WAKE,  SUPEREQUILIBRIUM  CALCULATIONS 
SUBROUTINE  WAKSC 
C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  COMPUTES  THE  AUXILIARY 

c turbulence  via  superlouilibrium  THEORY 

C 

C l)  FOR  THE  RIGHT  HaNU  SIDE  OF  POISSONS  EQUATION 
C 

C 2>  FOR  OU/OY,  OU/u Z,  DRHO/DY  AND  ORHO/qZ 
C 

C 3>  FOR  TURBULENCE  OUTPUT  TO  THE  LINE  PRINTER 
C 

DIMENSION  TURB(lC) , I UT  AL l 560 ) 

♦COPY  (LMWAK) 

etui valence  (Tukb(D « vmaT(1) > , (0Ivt,gmiii) ) 

EQUIVALENCE  (FF,GM(13) ) , <FT,GM(14) ) 

c 

louo  FORMAT  * //36H  S jPEKEUuJLIBRIuM  ERROR  HAS  OCCURRED » 2l 5 , 3Ei5 . 5 ) 
loci  FOh  MaT ( //42H  SuPEREUUlLIbRIUM  CORRECTION  TO  ZERO  NOISE, m, 

1 14H  TIMES  AT  STEP, 14) 

1002  F0RMaT(//16X,7hTURB  KE»8X,7HMEaN  KE , 8X , 7HKH0^RH0 ,6X , 

1 UhRHO^RHO  OUT  ,4X,11HMEAN  KE  OUT  « 8X  « 4HAREA  , 

<;  9X.7HPSI  MAX,faX,llHQISSlPATlON/10X,8El5.S) 

1003  FORMAT (1H0«15X,7HF  TOTAL, 9X,5HF  MAX,BX,8HF  CHANGE, 7X, 

1 10HDIVG  ERROR, 7X,4HLIFT/10X,5E15.S) 

1004  F0RMAT(//37h  TURBULENCE  CORRECTION  FOR  MAX  BOUNDS, 21b, 

1 AH  AT  STEP, 14) 

C 

FPRZ(aROM,ARG,ARGH)sFZM#(ARGM+FZS^ARG-FZR+aRGP) 
FPRY(ARGM,ARG,ARGH)=FYMM  ahgm+fys#arg-Fyr*argp ) 
SPRZ(ARGM«ARG,ARGR)=SZ!'HARGM+SZ^ARG  + SZP+AKGP 
SPRY (ARGM,ARG,ARGP)=SYM*ARGM+SY^ARGvSYP^ARGP 

c 

C ZERO  PERTINENT  VARIABLES 

c 

IOLAY r 7 

IF  (LSTFL)  a00,?0,50 
SO  MOOO=1 
LFL=U 
LW=0 
LWW-0 
LZFAFrO 

PLaNE=1,0/FLOAT(LYLFF+2)/F^0aT(LZLFF+2) 

DEPSN=0.0 

AREA=0.0 

DIVT=0.0 

CVSrO.O 

XFI=0.0 

FF=0,0 

FT  = 0,0 

CALL  SFVFH0.0,TURBX,10) 

CALL  SFVMVtROWG, TOTAL, 560) 

SCMX=4.C^YSCAL^YS0aL^ZSCAL/( YSCAL^YSCAL+ZSCAL^ZSCAL) 

ZEROV (6 )rSCMX 
FMAXV(6)=SCMX 
C 

C INITIALIZE  FOR  PASSING  THROUGH  PROFILE 
C 


i 
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IK0NA=2 
1K0WR=4 

IF  (LSTFL)  96.96,95 

95  IK0MR=2 

96  NRSTX=l\lRST*LZLFF 
NRENDsNRST+jENu-1 
NPOS^NRST-1 
CALL  WAKRR(NKST-l.ZV) 

CALL  WAKRR(NKST,ZFV) 

DO  2b0  NR=NRSTX,NKEND 

C 

r.  READ  three  surrounding  hows  and  test  for  selected  z in  DOMAIN 

c 

IF  (NP-NRST)  104.10A.103 
103  NHOS  = (\)R 

CALL  S>  VMv(Z.ZM,N*WZF) 

CALL  WAKMR (2.1) 

CALL  SFVMV(ZP,Z,NWWZF) 

Call  rakmR(3,2) 

CALL  WAKRK( NK+1 ,ZFV  ) 

GO  TO  105 
10H  ZM=Z+Z-Zp 
105  UZM=Z-ZM 
OZP=ZP-Z 
OZT=ZP-ZM 
CALL  WAKDG(Z) 

C 

c step  through  all  y points  computing  auxiliary  quantities 
c 

y=ygldv<  iyps-d 

YPsYOLUV( IYFS) 

CALL  WAKMP(iROwR,NR+l,IYPS-l,XPYV,l) 

CALL  WAKMP( IR0WR,NR, IYPS-1.XV,1 ) 

IF  (NR-NRST)  108,109,109 

108  CALL  'wARRF  ( XPYvfXI*iYV  ,2) 

GO  TO  110 

109  CALL  WAKMP(IRO^R,NR-l,IYPS-l,XMYV,l) 

HO  IYPSX  = IYPS+LYLFF 

CALL  WAK'IPJ  1R0^T,NR,  1YPS-1.TV.1) 

CALL  WAKMP( IR0wR,NR+l, IYPS.XPPV.1) 

CALL  WAKMP( IR0„R,NR, IYHS.aZPV.I) 

IF  (NR-NRST)  113,114,114 

113  CALL  WAKRF(XPPV,Xf'PV,2) 

GO  TO  lib 

114  CALL  WAKMP( IRO*R , NR- 1 , I YPS , XMP V » 1 ) 

115  00  200  JY=IYPSX, IYPL 

ic  (UY-IYPS)  116,117,117 

116  CALL  WAKRF(XPPW,XRMV,1) 

CALL  LAKRF(XZPV,XZMV,1) 

CALL  WAKRF(XMPV,Xr-MV,l) 

YM=Y+Y’YP 

GO  TO  130 

117  CALL  WAKMP(IR0wT,NR,UY,TV,1) 

CALL  SFVMV(XPYV ,XPMV,NMOVE) 

CALL  SFVMV(XV,XZMV,NMOVE) 

CALL  SFVMV(XMYV,Xf’MV,NMOVE) 

YMrY 

Y = YP 

YP=YOLOV( JY+l) 


L 


f 


T 


, 


0 

u 

ii 

y 


Li 


h 

l i 

Li 


? ' 
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CALL  WA*MP( IR0*R,NR+1, JY+1,XPPV, 1) 
CALL  WAKMP( TROwR,NRt JY+ltXZPVtl) 

IF  (NR-NRST)  125,128, 128 
12b  CALL  WAKRF(XPPV,XMPV,2) 

GO  TO  130 

128  CALL  WAKMP(lROwR,NR-l,dY+l,XMPV*l) 
130  DYM=Y-YM 
UYP=YP-Y 
DYT=Yp-YM 
C 

C COMPUTE  SPACING  FACTORS 
C 

FZM=-DZP/UZM/UZT 
F^P=OZM/DZp/UZ1 
FZ=-FzM-FZP 
TEMrOZiVDZP 
FZRrTEll+TLM 
FZS=FZR-1.0 
SZM=2, 0/OZM/QZT 
SZt'=2.0/DZP/UZT 
S2=-SZM-SZP 
C 

fym=-dyp/uym/oyt 

FTprDYM/UYP/UYT 
FY=-FYM-FYP 
TEM=UYM/D YP 
rYK  = TEivl*TtM 
FYS=FYR-1.0 
SYM=2. O/OYM/OYT 
SYP=2.0/DYP/OYT 
SYs-SYM-SYP 
C 

C check  APPROPRIATE  DENOMINATOR  values 

c 

IF  ( X V ( 1 ) ) 135,141,141 

13b  XV(1)=0.0 

LPL=LFL+1 

141  Q=SQRT(XV( 1)  ) 

FKZ=FpRZ(XMYV(2) , X V t 2 ) ,XPYV<2>  ) 

ftem=frz+ukdz 
CALL  WAKcHXVtFTEh) 

IF  (XV(6) -SCALM)  1410,1411,1411 

1410  XV ( 6 ) =3C ALM 

1411  IF  (XV(6)-SCMxi  1R2U, 1420, 1412 

1412  XV ( 6 ) =SCMX 

1420  IF  (LAMIN)  142,143,142 

142  Q=0.0 

GO  TO  150 

143  BQL=BBETa*Q/SCalE 
ci=a*bol*q/scale-g*ftem 
C2=BBS*G*FTEM-Cl 

IF  (C1*C2)  150,14^,150 

144  IF  (0)  145,150,145 

145  WRITE  ( NOUT , 1000  ) NR , UY , XV ( 1 ) , C 1 , C2 
LTRNFrl 

RETURN 

C 

C COMPUTE  DERIVATIVES  AND  RIGHT  HAND  SIOE 
C 
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150  FRY  = FpRY  <XZMV<2) ,XV<2) ,XZPV<2> ) 

FUZ  = FPRZ(XMYV(o) » X V ( 3 ) , XP Y V ( 3 ) ) 

FUY=FPRY ( XZM V ( 3 ) , XV ( 3 ) , XZP V ( 3 ) ) 

TV ( 1 ) =FUY 
TV(2)=FUZ 
T V ( 3 ) =FR Y 
TV ( 4 ) =FRZ 
IF  ( Q ) 153,153,152 

152  TtM=FuY*FUY*SCALE/G/BUL-BBETA*G*FTEM/C2 

1 -FUZ*FUZ*8BETA*(G*t-TEM*(l»0-A/B/S)+A*A*XV(l)  /SC AlE/BC ALE  ) /C1/L2 
CTEM  = CW  + (B-CVV*TLM)  /(HBETA  + TEM  ) 

CALL  NAKTC(CTEM,1.0,LVV) 

153  IF  (NSTAT-l)  1540  * 1335, 1540 
1555  T V ( 5 ) = 0 • 0 

TV ( 6 ) =0  • 0 
FVZ=0.0 
FwY  = 0 . 0 
GO  TO  160 
1540  TEMF=TV<5) 

CALL  WAKLL(TEMK,2) 

IF  ( 0 ) 154,154-155 

154  CTEM=0.0 
TEM=U.O 
TEMC=0.0 
GO  TO  156 

155  C5=CTEM* (OROZ+b*BBS*FKY*PRY/Cl) 

TEM=2 . 0*cTEM*G*FRY*PPRZ ( FPKY ( XMMV ( 1 ) , XMYV ( 1 ) , XMP V ( 1 ) ) , 

1 FPRY(XZMV(1) ,XV(1) ,XZHV(1) ) , FPR Y ( XpMV ( 1 ) , XpY V < 1 ) , Xppv ( 1 ) ) )/Cl 
TEMC  = CTEM-2.0*G*<L3*CTEM*F,<Z)/C2 
CALL  U/AKTC(TEMC,1.0-CTEM,LWW) 

TEMC=TEMC*SPRZ(XMYV<1) , XV U ) , XPY V ( 1 ) ) 

15b  FVY=FPRY(XZMV(4) ,XV(4) , XZpV<4) ) 

FVZ  = FpRZ(XMYV(‘+)  ,XV(4)  ,XPYV<4)  ) 

FWY  = FpRY (XZMV(5) ,XV(5) ,XZPV(5)  ) 

FW2=FPRZ(XMYV<5) , X V ( 5 ) ,XHYV<5>  ) 

TEMA=2.0*FVZ*F«Y 
TLMB  = 2»  0*F VY*F*Z 

T V ( 5 ) = -TEMC-CTt.M*6PKY  ( XZi'lV  ( 1 ) , XV  ( 1 ) ,XZPV(  1)  ) -TEM-TEMA+ I EM8-G* FRZ 
DXX=DX 

IF  <UXX-OIVP*OXMAX)  157,156,158 
157  OXX=UI VP*DXMAX 
15b  DXX=DIVF*0XX 

TV(5)=TV (5)+U*(FVY+PWZ)/0XX-TEMK 

160  CALL  WAKMP( lR0rtT,NR,JY,TV,2) 

C 

C WRITE  COMPLETED  SOLUTIUN  FOR  NEXT  INTEGRATION  STEP 
C 

CALL  WAKMP< IROwA,NR,JY,XV,2) 

C 

C COMPUTE  RICHARDSON  SCA*-E  LENGTH  COMPARISON 
C 

IF  (YOUT-YM)  170, 161 » 1605 
1605  IF  (YOUT-Y)  161,170,170 

161  IF  (LZFAF)  170,162,170 

162  IF  <4.0*xV(l)-FMAXV(l) ) 163,163,170 

163  CWS=ABS(FTEM) 

LZ  F = 1 

C 

C COMPUTE  SUPEREQUlLlfiRIUM  TURbULENCE  VALUES 
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c 

170 

171 

C 

c 

c 

17* 


174 

174 

175 
17b 

176 
C 

C 

C 

180 

181 

182 

186 

18** 

185 


18b 

167 

186 

200 


205 

250 

C 

C 

C 


IF  <Q)  171.171*172 

Call  sfvfuo.o.turs.io) 

GO  TO  176 

OKbEK  - u V UW  \ln  UR  VK  Lr\  RR  UO  VV  W* 

TURB(9>=CTEf-i*Xv(i) 

TUKB<3)=TUKtj(9)*G*FHY/Cl 

TUKB(5)=-daL*TURB<3)/G 

TUKB(1)=-«TUR6(9)*FUY*TURB(3)*FUZ)/BQL 
TURb(6)  = (B0L*( TURB(9)*FTEM+TURB(3)*FRY  ) 

1 -TURB(b)*BH£TA*G*FKY/B/S)/C2 
TUKB(10)=TURB(9)-2.0*G*TURb(6)/BQL 
CALL  wAKtC(TOKB(10) .XVtl)-TURB<9> .1 ) 
TURB(B>=XV(l)-TURB(9)-TURB(10) 

TURB(7)=-SCALE*(TURB(5>*FKY+TURB(6)*FTEM)/4l/B/S 
TURb<4)=(FTEM*<  TUKB*6)*FUY  + TURB(10)*FUZ) 

1 -btlL*  < TURB  ( 1 ) *F  RY  + I URB ( 5 ) *FUY+TURB  ( 6 ) *FU* ) )/Cl 
TURB(2)=-(G*TURB«4)+TUKB(3)*FUY+TURB(10)*Fu*)/BUL 
TURB(3)=TURB(3>-U*SCALE*<FwY+FVZ)/3.0 
00  174  1=1.10 

IF  (ABS(TUKB(I))-AB5(TURBX(I) ) ) 174,174.17b 
TURSX ( I ) =TURB ( i ) 

CONT IN  JE 

IF  (CTEM-CVS)  i7b  * 17b , 175 
CVS=CTEM 

IF  (LSOLF)  178.160,178 

CALL  L'AKMP(  IR0*G,NR»  JY,TURb,2) 

COMPUTATION  OF  INTEGRALS 


SUMF=PLANL*L  'T*UZ  I 


IF  (NR-NRST)  x 3 1 , 1 82 , 162 
SUMF=0 • 5*SUhF 


IF  (JY-NYPS)  163,184. 184 
SUMF=0.5*SUMF 

DEPSN=OEPSN+U*XV ( 1 ) *SUMF/SC ALE 

IF  (NSTAT-1)  135.200,185 

XF I=XFI*TV<5)*SUMF 

arla=aRla+sumf 

TEM=FvY+FRZ 

oivt=divt+tlm*tem 

IF  (ABS(TV(S) )”AB5(FP)  ) 167,187,186 
FF=TV ( 5 ) 

TEM=ABS»TVl5)-TEMF ) 

IF  (TEM-FT)  200,200.188 


FTsTEM 

CONTINUE 


CALL  WAKmRI 
IF  < LSOLF) 
CALL  RmKMG ( 
CONTINUE 
CALL  SFVMV ( 


NR.ZV) 
205.250,205 
NR.ZV, 2) 

TOTAL, RORG, 560 


) 


COMPUTE  ANU  PRINT  PERTINENT  INTEGRALS 


EPSSV(2)=SQRT(RiS*FMAXV(l)/G/CRS)/2.0 
IF  (LW+LRW)  252.252.251 
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WRITE  ( NOUT  « 1004 ) LV V , LWW , NPTSN 
IF  (LFE)  260  • 260  • 25& 

WRITE  ( NOUT  1 1001 ) LFL, NPTSN 
OEPST=OEPST  + 0.b*DXSAV*(DEP&l4>DEPSN)*B 
OEPSIrDEPSN 

WRITE  (NOUT, 1002)  (6M( I ) , 1=1 , 7 ) ,OEPST 
IF  (NSTAT-1)  270  . 300 ,270 
IF  <FKAXV(4)*FhAXV(5) ) 271.272*271 
DIVT=DIVT*YSCAL*2SCAL/(AREA*FMAXV(4)*FMAXV(5) ) 
WRITE  (NOUT, 1003)  XF I , FF , F I , DI VT , GM ( fi ) 

RETURN 

END 

ID  0105  OB  aDDR  4AE0  DB  CNT  0264 
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WAKSE.S(0105) 

**WAKSE  - STRATIFIED  SUBMAKINE  WAKE,  SET  EPSILON  CRITERION 
SUBROUTINE  WAKSE(tPS,EPSV) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  TRANSFERS  EPS  TO  EPSV 
C 

OIMENSION  EPSV(b) 

(CMWAK) 


♦ COPY 
C 


lOO 


DO  110  I = 1 » N V A K 
TEM=1.0L10 
VWTF=VWTFV( I ) 

IF  (VWTF)  10Q.108,1U0 
VMAX=GMAXV(I ) 

VSCA=VSCAV( I) 

IF  (VFAX-VSLA*EPSM  101,10*, 102 
VMAX=VSCA*EPSN 
IF  (EPS)  1031101+,  lo4* 

tem=-eps*vsca 

GO  TO  106 

tem=eps*vmax 

TEM=TEM/VWTF 
EPSV ( I ) =T  lM 
CONTINUE 
RlTURN 
END 

cart  id  0105  ob  aodr  2730  ob  cnt  0032 


101 

13* 

lUi 


10*+ 

10b 

lOP 

110 
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WAKSR.S(OlOb) 

♦♦WAKSR  - STRATIFIEO  SUBMARINE  WAKE,  SPECIAL  RESTART  profiles 
SUBROUTINE  WAKSR  (I- ILIO) 

c 

C THIS  SUBROUTINE  IN  Th£  WAKE  PROGRAM  LOCATES  A SET  OF 
C SPECIAL  RESTART  PROFILES  ANO  INITIALIZES  THE  WORKING  FILE 
C 

DIMENSION  F I L I D ( 3 ) , RdUF (6,40 ) , IBUF ( 40 ) , I VEC (3) 

♦COPY  ( CMWAK ) 

EQUIVALENCE  ( R sUF (1,1)  ,R0WG( 1 ,1 ) ) , ( IBUF (2) ,ROWG( 1 ,16) ) 

C 

OATA  lVEc/3,5,b/ 

C 

lOOu  FORMAT(//25H  IMPROPER  SPECIAL  RESTART, 415) 

C 

C locate  proper  X relord 
C 

NR  = 2 

K=1VEC<LZFAF) 

2l  0 NP=NR-1 

CALL  >3B)-0R(FlLl0,NP,2,ZA> 

NY=IVECA 

NZ-IVECB 

NR=NP+NY+NZ+k*nY*NZ 
CALL  Pdl-OR(FILIO,NRtl,XP) 

IF  (XP)  215,212,212 
212  IF  (XP-X)  210,^10,215 
2l5  X=ZA 

IF  (LZLFF)  2151,2152,2152 

2151  NRST=2 

GO  TO  2153 

2152  NRST=(MXRZ-ImZ  ) /2  + 2 

2153  JEND=NZ-2 

TF  (l-YLFF)  2154,2155,2155 

2154  NYPS=2 

GO  TO  2156 

2155  NYPS=(MXKY-lMY)/2  + 2 

2156  NYPE=NY+NYPS-3 

IF  (JENO)  218,218,2160 
2160  IF  (NRST-2)  218,2165,2165 
2165  IF  (NYPS-NYPE)  2l70,Z18,2l6 
2l70  IF  (NYPS-2)  216,220,220 
218  WRITE  ( NOOT ,1000)  N YPS , NYPE , NRST , JENO 
CALL  EXIT 

22 0 CALL  P8FOR(FILiO,NP,NY,YOLDV(NYPS-l) ) 

CALL  PBf-QR  (F  ILI  J ,NP,NZ,ZOLDV  (NRST-1 ) ) 

C 

C BUILD  WORKING  FILL  FROM  SPECIAL  RESTART  PROFILES 
C 

CALL  SFVFL(0,0,TV,NWVEC) 

00  250  JZ=1,NZ 
00  2200  JY=l,4o 
IBUF ( JY)=0 

CALL  SFVMV(ZERUV,RBUF(1,JY) ,NVAR) 

RBUF (6, JY)=ZEROV( 7 ) 

2200  CONTINUE 

00  230  I VAR-1 , NwvEC 
JVAR=I VAR 
KVAR=IVAR 


T 


* 


2202 

2203 

2204 
2206 
2207 
2206 
221 


2214 

221b 

222 

22*+ 

22b 

230 


CaKT 


GO  TO  (221, 221, 2202, 2204,2204, 2203), IVAR 

IF  (IMSTAT-2)  2207,230,2207 

KVAR=KVAR+1 

:c  (NSTAT-l)  2206,230,2206 
IF  (LZFAF-2)  230,2208,2208 
IF  (LZFAF-2)  221,230,221 
JVAR=jVAR+LZFAF-3 
NR=NP+ ( JVAR- 1 ) *I\IY*N2 
00  225  UY=1,NY,10 
I=JY/10+1 
J=10* ( 1-1 ) 

NK,'=MMIIM(10,NY-J) 

NRXX=NR+NRY* ( J2-1  )+J*NZ 

call  pbforifiliu.nrxx.nry, turbx) 

DO  224  1=1, NRY 
U=U  ♦ 1 

RBUF( IVAr,J)=TURBX( 1 » 

TEM  = ABS(Rt)UF(  IVAR,  J)  ) 

IF  (ABS(RBUF(IVAR,J)-ZEROV(KVAR) )-1.0E-04)  2214,2214,2215 

IBUF( J)=IBUF( J)+l 

IF  (TEM-GMAXV(RVAR) ) 224,224,222 

Gi“IAXV(KVAR)  = TEh 

CONTINUE 

CONTINUE 

CONTINUL 

NK=NRSr+jZ-2 

NPOS=MK 

Z=ZOLDV(NR) 

IYPirn 

IYPE=0 

DO  240  NRY=1,NY 
jY=NYPS+NRY-2 
IF  (IYPB)  231,231,233 
IF  (IBUF(NRY)-K)  232,240,240 
IYPS=FiiwlAX(2,JY  ) 

IF  ( IBUF(NRY)-K)  234,241,241 
IF  (NSTAT-l)  235  ,236,235 
TV(6)=RBUF(6,NKY) 

CALL  WAKMP(IR0WT,NR,JY,TV,2) 

RBUF(6,NRY)=SCALE 

CALL  WAKMP(  IROwR , NR , UY , RdUF ( 1 , NRY ) ,2) 

CONTINUE 

IF  (1YPS)  242,242,216 

I YPE=jY-l 

CALL  WAK^R(NR.ZV) 

CONTINUE 

RETURN 

ENU 

10  0105  OB  AOOR  57A0  DB  LNT  0004 
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WA*SY.S(0105> 

**waksy  - stratified  submarine  wake,  step  solution  in  y 
SUBROUTINE  waksy 
c 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  MAKES  A STEP  OF  DELTA  X/2 
C WITH  Z DERIVATIVES  EVALUATED  AT  THE  PRESENT  X POINT 
C 

♦COPY  (CMWAK) 

c 

C INITIALIZE  COMPUTATION 

c 

I olay =5 

GO  TO  ( 10  «20 ) « NSS 
10  IROWA=3 

IK0WR=2 
MUOD=-l 
LPKRF  = 0 
GO  TO  30 
20  IR0WA=4 

IR0WR=3 
MUUD=0 

IF  (LPKRF)  180,30,180 
30  NK=NRST-1 

npos=nr 

CALL  WAKRR ( NR , Z V ) 

CALL  WAKRR(|\iR  + l,ZPV) 

JENDXrJENu 

IF  (LZLFP)  100,102, 102 
loo  JENDX= JEND+ 1 

ZM=Z+Z-ZP 
I YPEMrlYPEP 
IYKSMrlYPSP 

102  DO  157  JZ=1,JENDX 
JZN= JZ+LZLFF 

IF  (JZN)  104,104,103 

103  NR=NR+1 
NPOS=NR 

CALL  SFVMV(Z,ZM,NWWZF) 

CALL  WAK MR ( 2 , 1 ) 

CALL  SFVMV(ZP,Z,NWWZP ) 

CALL  WAKMR (3,2) 

CALL  WAKRR(NR+1,ZPV) 

104  IYPSNzIYPS 

iypen.iype 

IBOT=lYPS+LYLFF 

lTOP=lYHr 

C 

c initialize  for  thl  y direction  and  sweep 
C 

10b  Y=YOLDV( IYPSN-1) 

YP=YOLDV( IYPSN) 

CALL  WAKMP(IRUwR,NR+l, IYPSN-1, XPYV,1) 

CALL  WAKMP( IROwT,NR+l , IYPSN-1 , TP YV,1) 

CALL  WAKMP ( IROwR, NR, IYPSN-1, XV, 1) 

CALL  WAKMP ( IROwT,NR, I YHSN-1 , T V , 1 ) 

IF  (JZN)  1061,1081,1082 
1081  CALL  WAKRF(XPYV,XMYV,2) 

CALL  WAKRF(TPYV,TMYV,4) 

GO  TO  1083 


. i 


i 
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1002 

1003 


1007 

1008 
1009 


CALL  WAKMP( IR0WR,NR-1«1YPSN-1 ,XMYV,1) 
CALL  WAKMP{ iROwT .NR-1 , I YPSN-1 , TMYV , 1 > 
CALL  WAKMP(  iKOft*«MR  + lt  I YPSN , XPP V , 1 ) 
CALL  WAK*'iP<  IROwR.NR.  1YPSN.XZPV  ,1 ) 

CALL  WAKMP(lROwT.NR.lYPSN.TZPV.l) 

IF  (JZN)  1087.1007,1008 
CALL  WAKRF  (XPPv  » XI"PV  ,2) 

GO  VC  1089 

CALL  WAKMPt  IR0*R,NR-1,  IYPSN.  XI^PV  . 1 ) 
JY=IYPSN-1 

IF  (LYLFF)  109,110,110 
CALL  |«A*RF(XPPV,XFMV,1) 

CALL  WAKRF(XZPV,X2MV.l ) 

CALL  wAKRF(XMPv,XP'iv|V,1) 

CALL  WAKRF(T2PV,T2MV,3) 

YM=Y+Y- YP 
Gu  TO  120 


UPWARD  pa^0 
) JY=JY+1 

CALL  SFVMV(XPYv,XPMV,NI10Vl) 

CALL  SFVMV(XV.XZMV,l\nOVt) 

CALL  SFVMVfXMYV.XfMV.NMCVE) 

CALL  SFVMV  ( TV  . TZMV,iMHOVfc:) 

YM=Y 

Y=YP 

YPsYOLOV ( JY+1) 

CALL  WAKMP<IROwT,NR+1.JY,TPYV,1) 

IF  (JZN)  112.112.113 
1 CALL  WAKRF<  TPYv.TP' YV.4) 

GO  TO  114 

5 CALL  WAKMP(lRO*T, NR-1, JY, TMYV, 1) 

» CALL  WAKMP( IK0WR,NR+1,JY+1,XPPV,1) 

CALL  WAKMP( I KOwR , NR ♦ JY+1 » XZPV , 1 ) 

CALL  WAKMP ( I ROw T , NR . UY+1 . TZP V , 1 ) 

IF  <JZN)  116.116.118 
» CALL  WAKRF(XPPV,XrPV,2) 

GO  TO  120 

i CALL  WAKMP(IRO*R,l\sR-l,JY+l.XMPV«l> 

calculate;  Matrix  coefficients  anii  invert  for  gamma  ano  av 

j CALL  WAKMY 

DO  125  Isl«NWVtC 
TEM=1.0/(YMAT(I)-XMAT(I)*GM(I) ) 

GM< I )=TEM*ZMAT(I ) 

AV(I)  = TEM*(0wEC(  I)-XMAT(I)*AVU)  ) 

3 CONTINUE 

CALL  WAKMP<lRO*A,NR.UY,AV,2> 

CALL  wAKMP( IR0wG.NR.UY.GM,2) 

IF  (JY-IYPEN)  110.13U.130 

UPPER  BOUNOARY  CONOITION 


CALL  WAK£C(I) 

IF  (I)  131.140.131 

IYPEN=JY 

GO  TO  150 


i 


II 


< 


A-82 


add  points  to  the  righi  side  where  needed 


IF  (JY  + 2-MXRT)  141,141,170 
IF  (JY-NYPE)  110,148*148 
YOLDV( JY+2)=YP+DFKMX*( yp-y ) 
GO  TO  110 


initialize  for  downwaru  pass 


1500 


1507 


I YPEX=IYREN-1 

iypsx=iypsn+lylff 

IF  (IYPSX-IYPEN)  1500,1515.1515 
DO  151  I Y= I YPSX , I YPEX 

jy=iypex+iypsx-iy 

CALL  WAKMP( lROWA,NR«UY *XV ,1 ) 
CALL  WAKMP( IROwG,NR«OY,GM,1 ) 

DO  1507  1=1 , NWVEC 
XV(I)=XV(I)-GM(I)*AV(I > 

CONTINUE 

CALL  SFVHV(XV, AV, NWVEC ) 

CALL  WAKMP( IROwa,NR«UY, AV,2) 
CONTINUE 

IF  (LYLFF)  156,152,152 


151  CONTINUE 

1515  IF  (LYLFF)  156,152,152 
C 

C LOfcEK  BOUNDARY  CONDI l ION 
C 

152  DO  154  Irl.NVAR 

IF  C ABS< AV( 1 » -ZERUV ( I ) ) -EPSSV ( I ) ) 
154  CONTINUE 


154,154,160 


RETURN  TO  MAINLINE 


IYPS=FMAX(2, IYPSN) 

I YPE= I YPEN 

NYPS=MMlN(N\PS,iYFS) 
NYPE=MMAX(NYPE,IYFE) 
CALL  WAKWR(NR.ZV) 
continue 


monotonicity  CHECK  TO  ZERO  VALUES 


1570 


1571 

1572 


1573 


NRSTX=NRST+LZLFF 
NRLNU=NKST+JENu-l 
IF  (LYLFF)  1585,1570,1570 
1 = 0 

ll=lzlff 

IYPS=l2LFF»(MXRY-NYPS)+MXRY 

DU  1590  JZ=NRSTX ,NRLNO 

IF  (LL)  l57l ,1574*1574 

IF  ( IYPSV( JZ)-lYPS)  1572,1580,1579 

1 = 1 

IYPE=I YPS-1 
I YPS= I YPSV ( JZ ) 

DO  1573  JY=I YPS , I YPE 
CALL  WAKWZ( JZ-1 , JY ) 

CONTINUE 

IYPSV( JZ-1)=IYPS 
GO  TO  1580 


mm 


1574  IF  (IYPSV(JZ)-NYPS)  1578,1578,1575 

1575  IF  ( IYPSV( JZ)-IYP5)  1579,1580,1576 

1576  1=1 
IYPL=lYPSV( JZ)-i 

00  1577  JY=I YPS , I YPE 
CALL  WAKWZ(JZ,JY) 

1577  CONTINUE 

I YPSV ( JZ)=IYPS 
GO  TO  1580 
1576  LL  = -1 

1579  I YPS=I YPSV ( JZ ) 

1580  CONTINUE 

IF  (I)  1585,1565,1570 
1555  1=0 

LL=LZLFF 

IYPE=LZLFF*(1-!\iYPL)+1 
DU  1595  jZ=NRSTX,NRENO 
IF  ( LL  > 1586,1569,1589 

1586  IF  ( IYPEV< JZJ-IYPL)  1594,1595,1587 

1587  I=i 

IYPS=IYFE+1 
IYPE=I YPeV ( JZ ) 

DO  15C8  JY=IYPS,IYPL 
CALL  bAKtdZ  ( jZ-1 , j Y ) 

1586  CONTINUE 

IYPEV (JZ-1 JrlYPE 
GO  TO  1595 

1569  IF  ( IYPLV< JZJ-NYPL)  1590 , 1593 , 1593 

1590  IF  ( IYPE\M  JZ)-1YPL)  15^1,1595,1594 

1591  1=1 

IYPS=lYPEV(JZ)+l 

DO  1592  JY=IYP5,IYPL 
CALL  LimK^Z  < JZ,UY) 

1592  CONTlNUt. 

IYPEV ( JZ ) =IYPE 
GO  TO  1595 

1593  LL=-1 

1594  IYPErlYPLV( JZ ) 

1595  CONTlNUt 

IF  (I)  180,180,1585 
C 

C AUD  POINTS  rO  THE  LEFT  SIDE  HHEKE  NEEUED 
C 

160  IF  (IYPSN-2)  170,170,161 

161  IYPSN=IYPSN-1 

IF  (IYPSN-NYPS)  162,106,106 

162  YOLOV ( IYPSN-1 ) =YOLQV ( IYPSN ) -OFRMX* ( YOLOV ( I YPSN+1 ) -YOLQV (IYPSN) ) 
GO  TO  106 

C 

C reuoltion  of  number  Of  POINTS  REQUIRED 
C 

170  LPKKF=1 
180  RETURN 

END 

cart  id  oi05  oe  adok  4U50  ob  cnt  uicc 
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W A*SZ • S ( 0 1 05 ) 

♦♦WAKSZ  - STRATIFIED  SUBMARINE  WAKE,  STEP  SOLUTION  IN  Z 
SUBROUTINE  waksz 
C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  MAKES  A STEP  OF  DELTA  X/2 
C WITH  Y DERIVATIVES  EVALUATED  AT  THt  PRESENT  X POINT 
C 

♦COPY  (CMwAK) 

C 

C INITIALIZE  COMPJTmTION 
C 

IQLAY=2 

GO  TO  (20,10) , NSS 
10  IROWA=4 

IKOwR=2 
MOOD=-l 
LPKRFrO 
GO  TO  40 
20  IROwArH 

IKOWR=4 
MOOU=0 

IF  iLPKRF)  180,30,100 
30  NRSTNrNRST 

JENUNrUENU 

100  NK=NRSTN-l 

npos=nR 

CALL  wAKRR(NR.ZV) 

ibot=i/ps+lylff 

ITOP=IYPe 

CALL  WAKRR(NR+l,ZHV) 

IBOTPrIYPSP+LYLFF 
ITOPP=IYPEP 
JZ  = 0 

IF  (LZLFF)  101,1015,101b 

101  ZM=Z+Z-Zp 

I YPEMrIYPEP 
IYPSMzI YpSP 
GO  TO  110 
1015  NR=NRSTN 

102  JZ=JZ+1 
NPOS=NR 

CALL  SFVMV(Z,ZM,NWWZF) 

CALL  WAKMR (2,1) 

CALL  SFVMV(ZP,Z,NWWZF) 

CALL  WAKMR (4,2) 

ibot=ibotp 

I T0P  = I T OPP 

IF  (JZ-JENDN)  104,104,106 
104  CALL  WAKRR(NR+1,ZPV) 

IBOTPrIYPSP+LYLFF 
I TOPPzIYPEP 
GO  TO  110 
10b  IYPSPzO 
IYPEPrO 

ZP=Z+DFRMX^(Z-4M) 

ZOLDV(NR+l)zZP 
1 YPSV ( NR+1 ) =0 
IYPEV(NR+1)=0 
IBOTP=0 


I 
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HI 

i. 

1115 

J] 

H2 

n 

1125 

II 

112b 

Li 

113 

i 

u 

1135 

114 

J Y 

u 

ii 

1141 

D 

1142 

ii 

1144 

i i 

• v 

1145 

*- 

1 1 4fa 

1147 

-- 

1140 

**  *> 

1149 

ITOPPrO 


INIT 1 ALIZL  FOR  THl  Y QIKECTION  AND  SWEEP 


Y=YULOV( IYPS-1) 

YP=YOLDV( IYPS) 

CALL  WAKMP(  lRUwR,NR  + l,IYPt>-l,XPYV*l ) 
CALL  WAKMP(  IROwT  ,I\iR  + 1,IYPS-1,TPYV,1) 
CALL  WAKMP(  I ROwR  , l\R  , I YPS-1 , X V , 1 ) 

CALL  WAKMP( IROwT ,KR, 1YPS-1, TV, 1) 

IF  (J2)  111,111,111b 
CALL  WAKRF(XPYV,Xr'YV,2) 

CaLL  wAKRF  ( TPYv  , Tl*  Y V ,4) 

GO  ro  112 

CALL  WAKMP(lROwR,KR“l,IYPS”l«XMYV*l) 
Call  WAMP(  I ROwT,KR-l,  IYP5-1,TMYV,1) 
lYPSXrlYPb+LYLFF 
IF  (J2)  112b, 1126, 1125 
CALL  WAKMP(IROwA,NR-l,IYPa’-l,AV,l) 
CALL  WAKMP ( lKO»b,NR-l, IYPS-1,GM,1) 
CALL  WAKMP ( IROi*R  , IViR  + l , I YP5 , XPPV * 1 ) 
CALL  WAKMP(  lROwR,l\R,  1YPS,XZPV,1) 

CALL  wAKMP(  IRO*T,IVR,IYPS,  TZPV,1) 

IF  IJZ)  113,115,1135 
CALL  WAKRF  ( xPPv  , X •wi P V , 2 J 
GO  TO  114 

CALL  WAKMP ( IROwR  , IVR-1 , 1YPS,XMPV,1) 
LZLRFrO 

IYPSNrMXRY+LYLFF* (MXKY-2) 

IYPENrl 

DO  14s  JY=IYPSX, I YPL 
IF  (JY-IYPS)  1141,1142,1142 
CALL  WAKRF ( XPPv , XPMV , 1 ) 

CALL  WAKRF (XZPV,X2MV,1) 

CALL  WAKRF(XMPV  ,XM«IV,1) 

CALL  WAKRFdZPv  , TZMV,^) 

YM=Y+Y-Vp 
GO  TO  120 

IF  (JZ)  1145,1145,1144 

CALL  WAKMP<iROwA,KR-1,JY,AV,1) 

CALL  wAKMP(1ROwG,KR-1,JY,GM,1) 

CALL  SFVMV(XPYV ,XPMV,NMOVE) 

CALL  SFVMV(XV,XZMV,IMP10VE;) 

CALL  SFVMV(XMYV  , Xl“  M V , NMO  V£. ) 

CALL  SPVMV(TV,TZMV,NMOVE) 

YM=Y 

Y=YP 

YP=YOLDV( JY+1) 

CALL  W AKMP ( IROwT  ,I\R+1,JY,TPYV,1) 

IF  (JZ)  1146,1146,1147 
CALL  WAKRF  (TPYV  ,T(“iYV,4  ) 

GO  TO  1146 

CALL  WAKMP(  IROwT  ,I\R-1  ,JY,TMYV,1 ) 

CALL  WAKMP  ( I ROwR  , l\R  + 1 , J Y + l » XPPV  , 1 ) 
CALL  WAKMP ( lROwR,l\R,JY  + l,XZPV,l) 

CALL  WAKMP  ( lROwT,l\R,  JY  + 1,  TZPV,1) 

IF  (JZ)  1149,1149,1150 
CALL  WAKRF (XPPV , XMPV ,2) 

GO  TO  120 
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11-0  CALL  WAKMPt IR0wR,NR-1, J Y+l , XMPV , 1 ) 

C 

C CALCULATE  MATRIX  COEFFICIENTS  AND  INVERT  FOR  GAMMA  AND  A V 
C 

120  CALL  WAKMZ 

DO  125  I = 1 , NW VtC 

TEM  = 1.0/(YMAT(l)-XMAT(I)*GM(in 

GM(1)=TLM*2MAT(I) 

AV(I)=TEM*(DVEC(I)-XMAT(I)*AV(I)  ) 

12b  CONTINUE 

CALL  W AKMP ( lRQwA,NR,JY  , A V , 2 ) 

CALL  wAKMP( IR0wG,NRi JYiGM,2) 

IF  (JZ-JENDN)  129«13U,13U 
C 

C UPPER  bOUNOaKY  CONDITION 


129 

IF  (JY-IYPSP)  1291,1292, 1292 

1291 

IF  (IYPSP-2)  145,145,130 

• 

1292 

IF  (JY-IYPEP)  145,145,130 

130 

CALL  LAKEC(I) 

IF  (I)  14b i 132 i 14b 

. » 

132 

IF  (JZ-UENDN)  134 , 133  ,133 

133 

lzlrf=i 

’ ? 

1*4 

iypsn=mmin(jy,iypsn) 

IYPSN=MMAX(2,IYPSN) 

» • 

IYPLN=MMaX( jY,IYPtN) 

IF  (JZ+2-MXRZ)  145,170,170 

'1 

C OUlPUl  CURRENT  ROw  TO  uisk 


c 

14b 

CONTINUE 

. i 

CALL  WAKwR’NR,ZV) 

CALL  WARMS ( NR , Z V , 2 ) 

IF  (JZ-JENDN)  147*146, 14b 

. I 

14b 

IF  (LZERF)  1465,150*1465 

I4b5 

IYPSPslYPSN 

r l 

GO  TO  1485 

C 

. i 

C AUD 

p 

i POINTS  TO  THE  LEFT  Slut  WHERE  NEEDED 

* 5 

L 

147 

IF  ( IYPSN-IYPSP)  1475,146,140 

, * 

1475 

IYPSP=IYPSN 

c 

C AUD  POINTS  TO  THE  RIGHT  SIDE  WHERE  NEEqEu 

c 

148  IF  (IYPEN-IYPEP)  149,149,1485 
148b  IYPEP=IYPEN 

149  NR=NR+1 
GO  TO  102 

C 

C INITIALIZE  FOR  DOWNWARD  PASS 
C 

150  NR=NR-1 
JLNUNrJZ 

I2END=JENON-LZLFF-l 
DO  154  IZ=1,IZEND 
NPOSrNR 

CALL  SFVMV(ZiZP*NHWZF) 

CALL  WAKMR(2*3) 


J 


I 


1 

u 

I! 


! 1 


CaLL  WAKRK(NKt^V) 

CALL  Vn AKMG  (iMR  « 2 V « 1 ) 

IdOT=lYPS+LYLFF 

IToP=IYP£ 

IYPSN=ivIXRY+LYLFF*(MXKY-2  ) 

IYPlN=1 

IYc.t\ID  = IYPLP-I  YHSP-LYLFF  + i 
DU  15Q3  IY  = 1 « I Yc.NL 
J Y = 1 YpEP- 1 Y + 1 

call  wamp<irowa,iyr+i,jy,av,1) 
call  lakmp<  iro*a,nk,  jy,x\z,u 

IF  (JY-1YPL)  1500,1500,1502 

1500  It-  tJY-lYPS)  1501 ,1506 , 150b 

1501  IF  (1YPS-2)  150b  , 1506 , 1502 

1502  DU  15o«5  I=1,NVAK 

If-  (ABS(AV(I|-^LRUV(I))-lPSSV(D)  1504,1505,150b 

1503  CONTINUE 
GO  TO  lbob 

1505  IYPSN=MMIN( jY, iYPGN) 

IYPS>N=MMaX(2»IYPSN) 

I YPEN-MMaX ( jY , 1YPLN) 

GO  TO  lbob 

1506  CALL  InAKMP(1R0w&,NR»JY,GT1»1) 

DO  1507  I = 1,|\|WVEC 

XV ( I ) =XV ( I ) -Gi“l  ( I ) *AV ( I ) 

1507  CONTINUE 

CALL  WAKMP(  IR0WA,NR, JY ,XV,2) 

1508  CO^lINUt 
C 

C ADD  POINTS  TO  THE  LEFT  SIjE  WHERE  NEEuED 
C 

IF  (IYPSN-IYPS)  lblb,152,ib2 
1515  IYPS=IYPSN 
C 

C AUu  POINTS  to  the  RIGHI  Slut  WHERE  NEEqEU 

c 

15c  IF  (IYPEN-IYPE)  153,153,1525 

1525  IYPt=lYPEN 

155  CALL  WAKwR(UR*Zv) 

NKrNR-1 

154  CONTINUE 
C 

C LOwEK  BOUNDARY  COnDIUON 
C 

IF  ( LZLFF ) 156 , i5H5» 1545 
154b  NR=NR+1 
L2EKFrO 

I YPSN=MXRY+LYLFF* (NxRT-2 ) 

IYPEN=1 

iyeno=iype-iyps-lylff+i 
DO  155  1Y=1«IYlMu 
JY=I YRL-IY+1 

CALL  LAKiv|P(iROwA,NR,JY  , AV  , 1 ) 

DO  1546  1 = 1 , NV  AR 

IF  (A8S(AV(I)-2ER0V(I)  )-EPSSV/(I)  ) 1 54b , 1 54b  , 1 547 
I54fe  CONTINUE 
GO  TO  155 

1547  IYPSN=MMIN( jY, 1YP5N) 

I YPbNrl'HlAX  (2,1  YPSN  ) 
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irPEN=MMAX< JT « 1TPEN) 

LZERFrl 

15&  CONTINUE 

IF  (LZEKF)  160  t 156  « 160 
C 

C RETURN  TO  MAINLINE 
C 

156  JENU=JEN0N 
NRSTrNASTN 
RETURN 
C 

C NEW  LOWER  ROW  REQUIRED 
C 

160  IF  (NRSTN-2)  170»A70»161 

161  NRSTN=NRSTN-1 

jenon=jendn+i 

Call  WAKRR(NR-1,ZPV) 

IBOT=o 

ITOP=0 

irpSM=irpsN 

I YPEMr I Y PEN 

CALL  WAKWR(NR-1,ZMV) 

ZOLUV (NRSTN-1)=2M-DF«MX*(Z-ZM) 
IYPSV(NKSTN-1)=0 
ITPEV(NRSTN-1)=0 
GO  TO  100 
C 

C REDUCTION  OF  NuMBER  UF  POINTS  REQUIRED 
C 

170  LPKRF=2 
180  RETURN 
END 

CART  ID  0105  DO  AODR  4F20  DU  CNT  0216 
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WAMC.S(0105> 

♦♦WAKTC  - STRATIFIED  SUBMARINE  WARE*  TURBULENCE-  CHECK  FOR  MIN/MAX  VALUES 
SUBROUTINE  WAKTC( TEHM, TMAXtLFL) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  CHECKS  THaT  VV  AND  WW  SATISFY 
C RATIONAL  BOUNDS  Oiv  ThEiK  BEHAVIOR  IN  AN  UNSTAdLE  DENSITY  GRADIENT 
C 

♦CUP Y (CMfcAK) 

C 

IF  (TEMM)  10»lu0»*0 
10  TLMM=0.0 

GO  TO  40 

20  IF  (TEMM-TMAX)  10U , 1U0  » 30 

3 0 TEMM=TMAX 

<+0  LFL  = LFL+1 


100  RLTURN 
ENU 


Cart  id  0105  db  auor  4220  db  cnt  0020 


L 


1- 
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wa*nwR.S(  0105) 

**WAKwR  - STRATIFIED  SUBMARINE  WAKE.  WRITE  A Z BUFFER  ROW 
SUBROUTINE  WAKwR(NRX.ZPOS) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  WRITES  A Z ROW  AND 
C FLUSHES  THt  ROw  BUFFER  AT  THE  NRX  POSITION  FOR  ALL  Y 
C 

DIMENSION  ZPOSU)  ,RLCB(2A) 

♦COPY  (LMWAK) 

EQUIVALENCE  ( KECb  < 1 ) « OVEC ( 1 ) ) 

C 

DATA  .JERRX/2HWK/ 

C 

C WRITE  ITPS  AND  IYPL  values  BmLK  from  zpos 

c 

IF  (NRX ) 300.300 ,10 
10  IF  (NRX-MXRZ)  20,20,300 

20  Call  SFVMV ( ZPOS , ZA , NWWZF ) 

ZOLUV (NRx )=ZA 
IYPSV(NRX)=IVECA 
IYPEV(NRX)=IVELB 
C 

C LOCATE  NRX  POSITION  1 N ROW  BUFFER  ARRAY 
C 

J=NRX-NPOS+t 
IF  (J)  S00 ,300 , 30 
3 0 IF  ( J-3 ) AO, AO, SOU 

Au  IF  (IVECA)  l50,15U,50 

50  I YPSX=IVECA+LYLFF 

I YPEXrlVECB 
C 

C WRITE  NRX  INFORMATION  FROM  BUFFER  ROW 
C 

NRXX=( NKX-1 )*MXRY+I YPSX 
DO  130  IrlYPSX, 1YPL* 

CALL  PBFDR( SLNID ,NRXX ,NWR, RECB) 

CALL  SFVMV (ROW3 (1,1 ,U),RLCB(1),NWVLC) 

IF  (MOOD)  105.110,115 

10b  CALL  SFVMV (R0Wb(7» 1 , U ) ,RLCB(7) , NMOVE ) 

GO  TO  120 

110  CALL  SFVMV(R0Wd(7. I »U) ,RLCb(13) ,NMOVL) 

IF  (I-IBOT)  112,111,111 

111  IF  ( I -I  TOP ) 120,120,112 

112  CALL  SFVMV(ZERuV,RLLB(7) ,NWVLC) 

GO  TO  120 

115  CALL  SFVMV(R0WU(7«I«U)  , RECB  ( 7 ) , NWVLC  ) 

CALL  SFVMV (K0WB(  13, 1 , J) , RECB ( 19 ) , NW /EC ) 

120  NRXX=NRXX-1 

CALL  P6FoW(SLN1G,MRXX,NWR ,RECB) 

130  CONTINUE 

150  return 

300  ULRR= jERRX 

RETURN 
END 
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WAKmZ , S ( 0105 ) 

**WaKwZ  - STRATIFIED  SUBMARINE  WAKE,  WRITE  A ZERO  POINT 
SUbROUTlNE  WAKwZ(NRZ,NKY) 

C 

C THIS  SUBROUTINE  IN  THE  WAKE  PROGRAM  WRITES 
C A ZERO  POINT  TO  THE  WORKING  FILE 
C 

DIMENSION  RLCB(24> 

♦CUPY  ( CMm AK ) 

EQUIVALENCE  (RlCB( 1 ) «DVEC< 1 ) ) 

C 

DATA  JERRX./2HWZ/ 

C 

If-  (NRZ)  100,100,10 
10  IF  (NRZ-MXRZ)  20,20,100 
20  NRXX=(t\lRZ-l)*MXRY  + NRT 

CALL  SFVFL(0.0,RECB,NWVEC ) 

DO  30  1=7 , NWR ,NwVtC 

CALL  SFVMV(ZLROV,RECB(I) ,NWVEC) 

30  CONTINUE 

CALL  PBFDW(SLN1D,NRXX,NWK,RECB) 
return 

100  JERRsJERRX 

return 

ENU 
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HAMN.SI  0101) 

♦IOCS (2501  READER, DISK) 

♦♦INITIALIZATION  PROGRAM  FOR  FULL  PLANE!  SWIRL  IN  WAKE!  PROGRAM 

define  file  u 32000,2, u,i\m) 

DIMENSION  YOLO V (40)  ,ZOLDV(40)  , VECUO) 

NINU=6 

READ  (NINU, 1000 ) NY 
1000  FORMAT(lA) 

READ  ( NINU ,1001)  ( YOLDV ( JY ) , JY=1 , NY ) 

10 Ui  FORMAT ( El2 • 4 ) 

R5=.4 
R V=1 , 5 
NZ=NY+NY-1 
DO  10  JY  = 1,ImY 
I YrNY-JY+i 
IZ=NZ-JY-H 

YOLOV(IY)=RS+YOLDV(lY) 

ZOLOV( JY)=-Y0LDV(1Y) 

ZOLOV ( IZ)=-ZOLuV ( JY  > 

10  CONTINUE 
NY  = NZ 

CALL  SFVMV(zOLOV, YOLUV.NY) 

RMAX  = ZOLDV/(NZ) 

NR=1 

RR=0.0 

WHITE  ( 1 * NR ) RR 

WRITE  ( 1 * NR ) NY ,NZ, (YOLDV(JY) ,JY  = 1, NY) , (ZOlDV(JZ) ,JZ“1 «NZ) 
DO  100  1=1. fa 
DO  90  JY=1,NY,10 
IY=MMIN( JY+9.NY ) 

DO  80  JZ=1,NZ 

CALL  SFVFH  0.0,V/EC,40  ) 

Z = ZOLDV(.jZ) 

IF  (JZ-1)  70,70,40 
30  IF  (JZ-NZ)  40,70,70 

HO  DO  faO  J=JY,IY 

Y=YOLDV( J) 

IF  (J-l)  60 , faO , 42 
42  IF  (J-NY)  4b , 60 , 60 

45  R=SuRT(Y*Y+Z*Z) 

IF  (R-RMaX)  46,60,60 

46  GO  TO  (51,60,52,54,54,60),! 

51  RR=(Z^Z+Y*Y)/RS/RS 

VEC( J)=.0l08/(1.0+RR/6,25)^*2 
GO  TO  60 

52  RK=Q.5^(Z*Z+Y^Y ) /RS/RS 
VEC(J)=0,080*(1.0-RR)*tXP(-RR) 

GO  TO  60 

54  RR=3,0*(Y*Y+Z*4) 

IF  (R)  60,60,5o 

55  IF  (R-RV)  555,60,60 

555  VEL=(1.0-EXP(-RR) ) ♦ ( KV -R ) ♦♦2/R/ 15 . 0 

IF  (1-4)  56,56,57 

56  VLC( J)=-VEL^Z/R 
GO  TO  60 

57  VEC( J)=VEL^Y/R 

60  CONTINUE 

70  WRITE  ( 1 * NR  ) (V/ECU')  ♦ J=JY,1Y) 

80  CONTINUE 
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